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1.  INTRODUCTION 


In  communications  systems  analysis,  one  deals  with  the  characterization 
of  data  for  which  the  underlying  process  occurs  naturally.  The  development 
of  natural  phenomena  is  usually  monitored  with  the  passage  of  time.  For 
example,  the  surface  temperature  of  the  ocean  is  measured  by  a ship  traveling 
in  a straight  line;  rainfall  density  is  monitored  as  a function  of  rate/hour; 
ionospheric  effects  in  communications  are  mapped  as  a function  of  time  of 
day.  Therefore,  in  analyzing  and  modeling  such  systems,  one  must  regard  data 
in  the  form  of  a time  series.  A time  series  is  a random  or  non-deterministic 
function,  x,  of  an  independent  variable,  t.  In  most  situations,  t will 
represent  time  or  some  other  physical  parameter  such  as  space.  The  charac- 
teristic feature  of  time  series  is  that  future  behavior  can  be  closely 
estimated,  but  cannot  be  predicted  exactly  as  would  be  the  case  for  purely 
deterministic  functions.  In  many  applications  of  analysis  and  modeling,  it 
is  convenient  to  assume  that  certain  physical  processes  can  be  described  by 
deterministic  functions.  However,  if  the  underlying  process  is  "stochastic", 
then  the  deterministic  point  of  view  may  give  misleading  information  with 
respect  to  the  independent  variables. 

The  field  of  statistical  analysis  as  applied  to  communications  is 
premised  in  the  fact  that  classical  functional  analysis  is  not  adequate  in 
dealing  with  a random  process.  A great  deal  can  be  deduced  from  those  random 
processes  that  are  stationary,  i.  e.,  a process  which  is  in  statistical 
equilibrium.  However,  for  those  physical  processes  which  are  non-stationary, 
analysis  may  be  untenable  and  the  "randomness"  or  non-stationarities  must  be 
dealt  with  in  order  to  bring  the  process  into  statistical  equilibrium  prior 
to  analysis.  If  this  randomness  is  not  properly  approached,  then  meaningful 
results  of  analysis  would  be  very  difficult  to  obtain.  One  example  of  a non- 
stationary  process  is  the  adverse  effects  on  High  Frequency  (HF)  radio  commu- 
nications caused  by  sudden  variations  in  the  ionosphere. 

The  question  of  characterizing  the  ionosphere  for  HF  radio  communica- 
tions has  been  approached  from  many  points  of  view,  j'l'J,  [2].  Primarily, 
information  gathered  over  twenty  or  so  years  is  analyzed  by  the  National 
Bureau  of  Standards,  £2j,  and  monthly  and  yearly  summaries  in  the  form  of 
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world  contours  are  published.  Ionospheric  forecasts  are  considered  to  be 
either  long-term,  short-term  £2],  or  real-time.  Long-term  forecasts  usually 
predict  undisturbed  monthly  median  conditions  at  a particular  hour  for  some 
specified  month.  They  may  be  prepared  to  cover  a long  period  of  time  in  the 
future,  i.e.,  one  year,  or  even  an  entire  solar  cycle  (22  years).  Long-term 
forecasts  are  most  useful  in  planning  and  management  of  the  HF  spectrum. 
Short-term  forecasts  usually  predict  ionospheric  conditions  in  the  near 
future.  They  are  prepared  by  modifying  long-term  predictions  by  using  values 
of  local  magnetic  activity  to  account  for  disturbances  caused  by  changes  in 
the  geomagnetic  field.  A true  real-time  prediction  scheme  would  probably 
require  that  forecasts  be  available  concurrently  with  ionosonde  observations. 

It  should  be  noted  that  the  ionosphere  is  composed  of  definable  layers 
of  differing  electron  density,  namely,  the  D,  E,  F1  and  F£  layers.  These  are 
subject  to  violent  and  random  changes  in  altitude  and  density  due  to  sunspot 
activity,  diurnal  (daily)  effects  of  the  sun's  radiation,  magnetic  storms, 
sudden  ionospheric  disturbances,  and  other  natural  phenomena.  These  changes 
may  either  occur  suddenly  with  little  warning,  or  they  may  take  place  in  a 
cyclical  manner  due  to  combination  of  these  layers  into  a single  F layer 
during  night-time  (diurna1  effects)  hours.  This  activity  of  the  layers 
directly  affects  the  reliability  of  HF  communications  to  the  extent  that  they 
may  cause  outages  for  extended  periods  of  time.  Due  the  random  nature  of 
this  phenomenon,  and  due  to  the  time  dependence  of  the  information,  a logical 
approach  to  forecasting  and  interpreting  the  results  with  respect  to  systems 
performance  seems  to  lie  within  the  realm  of  non-stationary  time-series 
modeling.  Therefore,  the  aim  of  this  report  is  to  develop  statistical  models 
to  forecast  in  near  real-time  and  to  characterize  the  underlying  stochastic 
process  of  short-path  oblique  incidence  (01)  and  vertical  incidence  (VI)  high 
frequency  (HF)  information  up  to  sixty  minutes  in  advance. 

In  Section  2,  a systematic  presentation  is  gi  ?en  of  the  essential  theory, 
philosophy,  and  modeling,  utilizing  the  time-series  methodology.  Section  3 
is  a presentation  of  the  modeling  and  analysis  of  500  Km  path  ionospheric  data 
acquired  between  Fort  Monmouth,  New  Jersey,  and  Fort  Drum,  New  fork.  A 
Spectral  Analysis  of  the  ionospheric  information  is  presented  in  Section  4. 
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2.  BASIC  CONCEPTS  IN  TIME  SERIES  AND  SPECTRAL  ANALYSIS 


In  this  section  the  basic  concepts  in  time  series  and  spectral  analysis 
are  presented.  We  explain  how  a time  series  can  be  thought  of  as  a reali- 
zation from  a stationary  stochastic  process  and,  hence,  be  described  by 
certain  statistical  functions.  The  conditions  that  insure  the  stability  of 
a linear  system  are  developed,  and  the  stationary  stochastic  models,  i.e., 
the  autoregressive,  the  moving  average,  and  the  mixed  autoregressive-moving 
average,  are  introduced.  We  also  develop  a "backward  filter"  whereby  these 
stationary  stochastic  models  can  be  used  to  describe  non-stationary  time 
series,  and  we  introduce  a procedural  approach  to  fit  the  models  to  non- 
stationary  time  series. 

With  regard  to  spectral  analysis,  a brief  and  basic  description  of  some 
of  the  concepts  that  will  be  utilized  in  the  analysis  of  climatological, 
ionospheric,  and  man/machine  interface  data,  will  be  given.  In  section  2.5 
we  shall  be  concerned  with  the  spectrum  in  general,  with  respect  to  the 
aforementioned  linear  stochastic  models.  In  section  2.6  we  shall  deal  with 
spectral  estimators  and  illustrate  the  manner  in  which  the  concept  of  the 
window  enters  the  scope  of  the  analysis.  The  cross  spectrum  is  defined  in 
section  2.7  along  with  its  properties  and  a brief  discussion  of  the  role  it 
plays  in  spectral  analysis.  Finally,  more  complete  treatment  of  the  window 
(types  and  properties)  will  be  given  in  sections  2.8,  2.9  and  2.10. 

2.1  TIME  SERIES  AND  STOCHASTIC  PROCESSES 

A time  series  can  be  thought  of  as  a sequence  of  highly  correlated 
successive  measurements  (serially  correlated)  representing  some  aspect  of  a 
physical  phenomena.  Each  of  the  measurements  is  associated  with  a moment  of 
time  and,  in  some  cases,  some  other  physical  parameter.  A time  series  can 
either  be  continuous  or  discrete,  depending  upon  whether  the  observations  are 
continuous  or  discrete.  In  this  study  we  shall  be  concerned  with  only 
finite  discrete  time  series  which  are  measured  at  equidistant  time  intervals, 
or  those  continuous  time  series  that  have  been  digitized  to  form  finite 

discrete  series.  We  shall  denote  such  a time  series  by  x... 
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Time  series  can  also  be  generally  classified  as  being  either  a 
deterministic  function  or  a non-determini  Stic  function  of  time.  A determi- 
nistic time  series  is  one  that  can  be  described  by  an  explicit  mathematical 
relationship-,  hence,  the  future  values  of  the  series  can  be  forecasted 
exactly.  Many  physical  phenomena  occurring  in  practice  produce  deterministic 
data,  such  as  the  increase  in  water  pressure  as  one  descends  into  the  oceans, 
or  the  path  of  a spaceship  to  the  moon.  However,  in  most  cases,  time  series 
occurring  in  practice  are  non-deterministic;  that  is,  they  exhibit  random  or 
fluctuating  properties.  Unlike  the  deterministic  time  series,  there  does  not 
exist  any  explicit  mathematical  relationship  with  which  to  forecast  exact 
values  in  the  future.  Hence,  for  those  time  series  which  are  random  in 
nature,  we  must  describe  them  in  terms  of  probability  statements  and  statis- 
tical averages  rather  than  by  explicit  relationships.  The  prime  area  of 
interest  in  this  study  will  be  in  non-deterministic  time  series;  hence,  it 
will  not  be  possible  for  us  to  precisely  forecast  future  values  of  the  time 
series. 

To  describe  these  non-deterministic  time  series,  we  use  the  concept  of 
a stochastic  process;  that  is,  we  consider  an  observed  non-deterministic  time 
series  as  a realization  of  a stochastic  process.  To  explain  the  relationship 
between  non-deterministic,  or  statistical,  time  series  ana  stochastic 
processes,  consider  the  following.  A given  time  series  (xt>  t s 1,  2,  ...,  n) 
representing  an  ordered  random  phenomenon  is  assumed  to  be  a single  sample 
from  a particular  generating  process  {Xt,  t = -<=°,  ...,  -1,  0,  1,  ...,«>}. 

This  collection,  or  ensemble,  of  all  possible  sample  time  series  which  the 
random  phenomenon  might  have  produced  and  its  associated  probability  distri- 
bution is  called  a stochastic  process.  Thus,  a given  time  series  from  a 
random  phenomenon  may  be  regarded  as  one  physical  realization  of  the  doubly 
infinite  set  of  functions  which  might  have  been  generated  by  the  stochastic 
process.  The  set  is  doubly  infinite  because  an  infinite  set  of  values  is 
possible  at  any  given  time  and  because  there  are  an  infinite  number  of  time 
points. 

A stochastic  process  is  said  to  be  strongly  stationary,  or  stationary 
in  the  strict  sense,  if  the  joint  probability  distribution  of  any  set  of 
observations  is  not  affected  by  the  shifting  of  all  times  of  observations 
forward  or  backward  by  any  integer  amount  k.  That  is,  if  the  joint  proba- 
bility density  function  associated  with  n observations  X ,X  , ...,X  made 

12  n 
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at  any  set  of  times  t-| , t2,  ...»  tp,  is  the  same  as  that  associated  with  n 
observations  X 1+k,  X2+k,  • ••,  Xn+k  made  at  times  t1+k,  t2+k,  ....  tn+k.  For 
a process  to  be  strictly  stationary,  it  is  necessary  for  the  entire  proba- 
bility structure  to  be  time  invariant. 

A stationary  stochastic  process  may  be  described  by  its  mean  y,  which 
can  be  estimated  by 


x 3 „ Ex,  (2.1.1) 

n t=l  r 

the  sample  mean  of  the  time  series,  and  by  the  variance,  a*,  of  the  stochastic 
process,  which  can  be  estimated  by 


n 

E (x 
t-1  Z 


x)2, 


(2.1.2) 


which  is  called  the  sample  variance  of  the  time  series.  As  mentioned 
earlier,  the  values  of  the  time  series  at  different  points  in  time  are 
serially  correlated.  This  correlation  is  of  great  importance  to  this  study; 
hence,  the  covariance  function  of  the  stochastic  process  is  of  great  impor- 
tance to  us.  The  covariance,  Yk>  between  xt  and  xt+k,  k intervals  of  time 
apart,  is  called  the  autooovariance  at  lag  k.  yk  = cov  (X^.,  X^+k)  can  be 
estimated  by 

, n-k 

cxX00  3 „ E (xt  - x)(xt+k  - x),  k = 0,  1 n-1  , (2.1.3) 

t— 1 


the  sample  autooovariance  function  of  the  time  series.  Of  equal  importance, 
or  possibly  greater,  is  the  autocorrelation  at  lag  k,  pk  = Yk/Y0>  which  acts 
like  a correlation  coefficient,  and  can  be  estimated  by 


rxx(k> 


cxx(k> 


0,  1, 


n-1 


(2.1.4) 


the  sample  autocorrelation  function  of  the  time  series.  Note  that  when 
k = 0,  rxx(0)  is  1.  Both  the  autocovariance  function  and  autocorrelation 
function  are  even  functions  because  of  the  stationary  assumption.  In 
practice  it  is  only  necessary  to  compute  the  sample  autocovariance  and 
autocorrelation  functions  for  lags  up  to  n/4. 
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One  of  the  most  important  assumptions  made  with  respect  to  a time  series 
is  that  the  corresponding  stochastic  process  is  stationary,  [3 j*  In  general, 
the  properties  of  a stochastic  process  are  time  dependent;  that  is,  the 
current  value  will  depend  only  on  the  time  which  has  elapsed  since  the 
process  began.  We  can  make  a simplifying  assumption  that  the  time  series 
corresponding  to  the  stochastic  process  has  reached  some  form  of  steady  state 
or  equilibrium,  in  the  sense  that  the  statistical  properties  of  the  time 
series  are  independent  of  absolute  time. 

2.2  STATIONARY  STOCHASTIC  PROCESSES 


Stationary  stochastic  processes  are  used  to  model  time  series  of  many 
practical  situations.  Consider  the  discrete  process  where  the  random 
variables  Z^,  t=l , 2,  ...,  n,  are  mutually  independent  and  are  normal  with 
mean  zero  and  variance  a|-- this  constitutes  the  simplest  form  of  a stationary 
stochastic  process.  For  this  process  the  autocovariance  function  is  zero  for 
all  lags  except  the  zeroeth.  Such  a sequence  of  random  variables  is  called 
a purely  random  process  or  white  noise. 

An  important  class  of  stochastic  processes  can  be  generated  by  the 
passing  of  a purely  random  process  through  a linear  system,  (see  Figure  2.1). 


WHITE  NOISE 


FIGURE  2.1: 


OUTPUT 


REPRESENTATION  OF  A LINEAR  SYSTEM 


LINEAR 

SYSTEM 


If  the  system  is  linear,  we  may  express  the  relationship  between  the  output 
process  Xt  and  the  input  process  Zt  as 

*t  - » ’ "k  h-k-  (2-2-D 
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is  called  the  weighting  function  or  the  impulse  response  function  of  the 
linear  system.  The  resulting  stochastic  process  derived  from  a purely  random 
process  using  (2.2.1)  is  called  a linear  process. 

The  variance  of  the  linear  process  is  given  by 

o o 

a5  = ai  l h*  (2.2.2) 

* L j«0  J 

The  convergence  of  the  above  series  ensures  that  the  linear  process  has 
finite  variance.  We  shall  now  show  an  equivalent  conditon,  [4],  [3],  which 
also  assures  the  variance  to  be  finite.  We  will  use  Z-transforms,  [5]»  to 
obtain  the  characteristic  equation  of  the  linear  process  and,  then,  by 
restricting  the  roots  of  the  characteristic  equation  to  a certain  region,  we 
will  be  in  a position  to  establish  the  conditions  necessary  to  ensure 
stationarity  and/or  invertibility  of  the  process. 

Consider  the  following  difference  equation: 


yr  = *lVl  + Vr-2  + * Vr-m  * boxr  * •"  * Vr-n  ■ (2'2'3) 

where  a, a . b . ....  b are  the  parameters,  and  x is  assumed  given 

1 m 0 n r 

for  all  values  of  r.  The  general  solution  to  (2.2.3)  is 

OO 

yr  ■ 2 h.x  k , (2.2.4) 

r k=Q  K r-x 

where  hk  is  as  defined  in  (2.2.1). 

The  x's  and  the  y's  in  equation  (2.2.3)  can  be  considered  as  being 
obtained  by  sampling  the  continuous  signals  x(t)  and  y(t)  at  the  moments  of 
time  t = rA,  (r-l)A,  ...,  (r-n)A,  and  t = rA,  (r-l)A,  ....  (r-m)A, 
respectively.  Of  course,  A represents  the  spacing  interval  used.  We  can 
now  rewrite  (2.2.3)  as 


y(t)  - a1y(t-A)  - a2y(t-2A)  - ...  - a(TJy(t-mA) 

* b0x(t)  + b^x(t-A)  + ...  + bnx(t-nA)  . 


(2.2.5) 
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The  Fourier  transform  of  (2.2.5)  is  given  by 


[1  - aie'J'2nfA 


V'^7^]  Y(f) 


= [bQ  + b.,e"J2TrfA  + ...  + bne"j27TfnA]  X(f)  . 
Solving  for  Y(f),  we  have 


(2.2.6) 


Y(f) 


[b  + b,e"^2lT^A  + ...  + b e'^27TfnA] 
- o n J 


[1  - aie-™A 


. a g-j^TrfmA-] 
m ■* 


(2.2.7) 


Hence,  the  frequency  response  function: 

H(f ) = 2 h.e"j2TrfkA  , 
k=0  K 

which  is  the  Fourier  transform  of  h,  , is  given  by 


H(f ) 


[bQ  + bl9'j27rfA  + ...  + bne'j2lTfnA] 

[1  - a,e'j27TfA-  ....  a e'J'27TfnA] 

m J 


(2.2.8) 


(2.2.9) 


The  concept  of  Z-transforms,  [5],  is  used  to  manipulate  the  frequency 
response  function  H(f).  Substituting 


z = ej2TTfA 

in  equation  (2.2.9),  one  obtains 

r-1 


(2.2.10) 


[b  + b^-1  + ...  + b Z'n] 
H(Z)  = —2 i S 


(2.2.11) 


[1  - a,Z  - ...  - a Z’m] 

I m J 

which  is  the  Z-transform  of  the  impulse  response  function  hk;  that  is, 


H(Z)  * Z h.Z_k  . 
k=0  K 


(2.2.12) 
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From  an  operational  point  of  view,  Z may  be  thought  of  as  a shift 
operator;  that  is. 


Hence,  the  linear  difference  equation  (2.2.3)  may  be  expressed  as 
(1  - a.|Z~ ' - ...  - amZ"m)yr  = (bQ  + ^Z"1  + ...  + bnZ_n)xr  . (2.2.13) 

Solving  for  y , we  have 


y 


r 


[1  - a,Z''  - • 


* V'"] 


- i,n 

m 


V » 


then  using  (2.2.11),  equation  (2.2.14)  can  be  written  as 


(2.2.14) 


yr  = H(Z)xp  . (2.2.15) 

H(Z)  is  called  the  transfer  function  of  the  discrete  linear  system. 
Expanding  H(Z)  in  powers  of  Z"1,  yields 

00  00 

yr  ‘ ^ xr  " ^ ^kxr-k  ’ (2.2.16) 

r k=0  k r k=0  k r k 

which  is  the  general  solution  of  (2.2.3). 

The  characteristic  equations  of  a linear  process  can  be  obtained  as 
follows: 

(1)  by  factoring  Z m out  of  the  denominator  of  (2.2.14),  substituting 
ip  for  Z and  equating  to  zero,  we  have 

V"  - a^"1'1  - ...  - am  ■ 0 . (2.2.17) 


A linear  process  is  said  to  be  stable  if  the  roots  r of  the 

i 2 m 

above  characteristic  equation  lie  Inside  the  unit  circle,  [3].*  If  this 

condition  holds,  the  process  is  said  to  satisfy  the  stationary  condition. 
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(2)  by  factoring  Z n out  of  the  numerator  of  (2.2.14),  we  have 


[bQZn  + b^"'1  + ...  + bn]  Z’n  . (2.2.18) 

Solving  equation  (2.2.15)  for  xr,  one  obtains 

xr  = H(Z)_1yr  • (2.2.19) 

Now  substituting  \p  = Z into  (2.2.18)  and  equating  to  zero,  we 
obtain  another  characteristic  equation  of  the  system, 

bQ«  + b^H"1  + ...  + bn  = 0 . (2.2.20) 

A linear  process  of  this  type  is  said  to  be  stable  if  the  roots  ^ , s-,..., 
?n  of  equation  (2.2.20)  lie  within  the  unit  circle,  [3].  If  this  condition 
is  satisfied,  the  process  is  said  to  be  invertible. 

Also,  note  that  the  invertibil ity  condition  is  independent  of  the 
stationarity  condition,  [4  ],  [3]. 

2.3  THE  STATIONARY  STOCHASTIC  MODELS 


Consider  the  special  case  of  the  linear  difference  equation  (2.2.3),  in 
which  the  first  p of  the  a's  are  non-zero  (p  <_ m),  bQ  = 1,  and  b^  * 0 i > 1. 
This  results  in  the  equation, 

yr  = Vr-l  + a2y r-2  + •••  + apyr-p  + xr  • t2-3-1) 

This  resulting  process  is  called  an  autoregressive  process  of  order  p.  We 
shall  write  the  finite  autoregressive  process  in  the  following  form: 

Xt  - u - o1(Xt_1  - u)  + «2(Xt_2  - u)  + ...  + ctp(Xt_p  - y)  + Zt.  (2.3.2) 

Note  that  y is  the  mean  of  the  process  and  Z^.  is  a purely  random  process. 

As  we  have  just  shown  in  the  preceding  section,  the  autoregressive  process 
can  be  thought  of  as  being  stable  by  satisfying  the  stationarity  condition 
that  the  roots  of  its  characteristic  equation  (2.2.17)  lie  inside  a unit 
circle. 
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Next,  consider  the  special  case  of  the  linear  difference  equation 
(2.2.3),  in  which  the  first  q of  the  b's  are  non-zero  (q  £ n)  with  bQ  = 1 , 
and  all  the  a's  are  zero.  Thus,  we  have  the  resulting  equation: 

»r,,rtblVl*-‘Vr-q'  . (2-3'3> 


This  process  is  called  a moving  average  process  of  order  q.  We  shall  write 
the  finite  moving  average  process  in  the  following  form: 

Xt  " u = Zt  ' SlZt-l  * ‘ 6qZt-q  ’ (2.3.4) 


Thus,  the  moving  average  process  is  said  to  be  stable  if  it  satisfies  the 
invertibil ity  condition;  that  is,  the  roots  of  its  characteristic  equation 
(2.2.20)  lie  inside  a unit  circle. 

We  now  introduce  a backward  shift  operator  which  will  be  very  useful  in 
manipulating  stationary  stochastic  models.  The  backward  shift  operator  B is 
defined  by 


(2.3.5) 


Using  the  backward  shift  operator,  it  will  be  shown  how  a finite  auto- 
regressive process  can  be  expressed  as  an  infinite  moving  average  process. 
Consider  a first  order  autoregressive  process  with  y * 0, 

*t  s alVl  * zt  * (2-3-6) 

or  using  the  backward  shift  operator  (2.3.6)  becomes 


(1  - a-|B)xt  = zt  . (2.3.7) 

Expressing  the  x's  in  terms  of  the  z's,  we  have 

x^  = (1  - a^B)*^  zt  . (2.3.8) 
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Expanding  (1  - c^B)'1  we  have 

*t  = 2t  + + a2zt-2  + ‘ (2.3.9) 

Equation  (2.3.9)  is  an  infinite  moving  average  process.  Likewise,  one  can 
show  that  the  moving  average  process  can  be  expressed  as  an  infinite  auto- 
regressive process. 

The  point  that  has  been  brought  out  in  the  above  paragraph  is  that  it 
may  be  necessary  to  include  parameters  from  both  the  autoregressive  and 
moving  average  processes  in  order  to  achieve  parsimony-- employing  the 
smallest  possible  number  of  parameters  for  adequate  representation,  [ 6 ]* 

Of  course,  this  is  due  to  the  fact  that  a moving  average  process  could  not 
be  parsimoniously  represented  using  an  autoregressive  process,  and 
conversely,  an  autoregressive  process  could  not  be  parsimoniously  represented 
using  a moving  average  process.  Hence,  we  will  have  the  stationary 
stochastic  model , 


(t  - u = ai(xt-l  " h)  + a2^Xt-2  " + •••  + ap^Xt-p  ~ 


* h - 3lZt-l  - - Vt-q  ' 


(2.3.10) 


Equation  (2.3.10)  is  called  the  mixed  autoregressive-moving  average  process 
of  order  (p,  q).  Writing  (2.3.10)  in  terms  of  the  backward  shift  operator  B, 
with  u zero,  we  have 


0 - o^B  - ...  - apBp)Xt  = (1  - - ...  - SqBq)Zt  , 

(1  - S,B  - ...  - S Bq)  (2.3.11) 

X = ! a i 

(1  - a,B  «...  - apBP)  t 

Hence,  the  mixed  autoregressive-moving  average  process  can  be  thought  of  as 
the  output  from  a linear  system,  whose  impulse  response  function  is  the 
ratio  of  two  polynomials,  when  the  input  is  white  noise  1^. 

As  previously  mentioned*  the  stationarity  and  invertibil ity  condition  are 
independent  of  each  other;  thus,  we  simply  have  to  satisfy  the  stationarity 
* Chapter  1 
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condition  for  the  autoregressive  process  and  the  invertibi 1 i ty  condition  for 
the  moving  average  process  included  in  the  mixed  process  for  the  process 
(2.3.10)  to  be  stable. 


2.4  PROCEDURAL  APPROACH  FOR  FITTING  FORECASTING  MODELS 
TO  NON-STATIONARY  STOCHASTIC  REALIZATIONS 

At  present  there  exist  techniques  to  analyze  stationary  time  series 
records.  However,  the  techniques  available  for  the  analysis  of  non- 
stationary time  series  are  inadequate  and  do  not  lend  themselves  to  mean- 
ingful interpretations  of  physical  problems.  It  is  possible,  however,  to 
adjust  non-stationary  time  series  so  as  to  be  able  to  apply  the  existing 
techniques  of  stationary  time  series  analysis  directly.  This  adjustment 
takes  the  form  of  applying  a proper  filter  to  the  observed  non-stationary 
time  series  to  filter  out  the  non-stationary  components . (See  section  3.3.1). 

In  this  section  we  shall  illustrate  the  procedure  to  identify  whether 
an  observed  series  exhibits  stationary  or  non-stationary  properties.  If  the 
time  series  is  non-stationary,  we  explain  how  it  can  be  filtered  and  intro- 
duce the  concept  of  a "backward  filter."  We  give  a procedural  approach  to 
determine  the  stochastic  model  which  gives  the  best  fit  to  the  observed 
series,  and  we  apply  diagnostic  checks  to  determine  the  goodness-of-fit. 
Finally,  we  discuss  how  the  fitted  stochastic  model  can  be  employed  in  fore- 
casting and  updating. 

2.4.1  Identification  and  Filtering 

In  a given  physical  situation,  one  will  have  available  a stochastic 
realization  x^ , x^,  ....  xn>  of  n observations.  The  first  concern  is  to 
identify  whether  the  series  xt  exhibits  stationary  or  non-stationary 
properties.  To  accomplish  this,  we  make  use  of  certain  statistical  methods 
in  conjunction  with  the  properties  of  stationary  time  series.  A stationary 
time  series  will  have  the  following  properties: 

i)  it  will  be  in  a steady  state  in  the  sense  that  it  is  in 
equilibrium  about  a constant  mean  level; 
ii)  it  will  contain  no  trend; 
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i i i ) its  sample  autocorrelation  function  will  dampen  out  rapidly. 
Hence,  we  will  first  plot  the  series  as  an  aid  in  exercising  some  judgment 
about  the  behavior  of  the  information.  Next,  we  shall  apply  various  non- 
parametric  trend  tests  to  the  data.  Finally,  and  of  greatest  importance,  we 
calculate  the  sample  autocorrelation  function  (2.1.4)  of  the  observed  data. 
From  these  data-analysis  tools,  we  will  have  sufficient  information  to 
identify  the  given  observed  series  as  stationary  or  non-stationary. 

If  we  identify  the  observed  series  as  exhibiting  non-stationary  proper- 
ties, we  need  to  find  a filter  that  will  remove  the  non-stationary  components. 
One  of  the  most  used  and  most  efficient  methods  of  removing  non-stationary 
components  from  a time  series  is  by  differencing.  A first-order  difference 
filter  is  defined  by 


yt  = xt  - *t_-|  . (2.4.1) 

where  xt  is  the  observed  non-stationary  series  and  yt  is  the  resulting 
filtered  series.  Similarly,  a second-order  difference  filter  is  defined  by 

wt  = xt  - 2xt_1  + xt_2  , (2.4.2) 


and  so  on. 

Since  we  will  be  almost  exclusively  dealing  with  non-stationary  time 
series,  and,  since  either  a first-order  or  second-order  difference  filter  is 
usually  sufficient  to  transform  most  practically  occurring  non-stationary 
series  into  stationary  ones,  [6],  the  presentation  shall  be  confined  to 
determining  the  degree  of  differencing  necessary  to  result  in  a stationary 
series.  Using  the  backward  shift  operator  B (2.3.5),  we  express  the  differ- 
ence filter  in  the  following  form: 

yt  * 0 - B)d  xt  . (2.4.3) 

We  must  determine  a suitable  value  for  d,  either  0,  1 , or  2;  zero  will 
correspond  to  the  fact  that  our  observed  information  is  stationary;  one  will 
correspond  to  the  fact  that  a first-order  difference  filter  is  necessary  to 
filter  the  observed  series,  and  so  on. 

The  procedure  to  identify  the  proper  value  of  d is  to  compute  the  first 
and  second  differences  of  the  observed  series  xt,  t 3 1 , 2,  . . . , n.  That  is, 
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Xj.  is  processed  through  a first-order  difference  filter, 

yt  = (1  - B)xt  = xt  - xt_1  , 

which  will  have  (n-1)  values  and  a second-order  difference  filter, 

wt  = (1  - B)  *t  = - 2xt_i  + xt_2  , 

which  will  have  (n-2)  values.  For  the  observed  xt  and  the  filtered  yt  and 
wt  series,  we  calculate  the  sample  autocorrelation  functions  and  conduct 
trend  tests. 

By  examining  the  sample  autocorrelation  function  and  the  result  of  the 
trend  test  for  the  separate  series,  one  should  be  able  to  infer  a suitable 
value  for  d,  specifically,  the  degree  of  differencing  necessary  to  induce  the 
sample  autocorrelation  function  to  dampen  out  fairly  rapidly  and  to  cause  the 
trend  test  to  be  insignificant.  This  will  yield  a stationary  series  with 
which  to  continue  the  analysis. 

It  must  be  noted  that  in  some  instances  the  first-order  and  the  second- 
order  difference  filter  may  fail  to  remove  the  non-stationary  components. 

When  this  occurs,  we  must  continue  to  search  for  a proper  filter  that  will 
leave  us  with  a stationary  series.  One  alternative  is  to  apply  a higher- 
order  difference  filter,  or  we  can  try  other  types  of  filters.  Jenkins  and 
Watts,  [3],  list  several  other  types  of  filters.  In  some  respects,  filtering 
a non-stationary  series  is  a trial  and  error  procedure  in  that  one  attempts 
to  transform  the  observed  series  into  a stationary  one  by  the  use  of  a 
mathematical  function.  Unfortunately,  this  requires  us  to  search  for  the 
proper  function  to  accomplish  this  purpose. 

From  the  examination  of  the  sample  autocorrelation  function  and  the 
sample  partial  autocorrelation  function,  which  is  defined  in  terms  of  the 
sample  autocorrelation  function  as 

i(U)  - rxx(l)  , (2.4.4) 
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and 


(2.4.5) 


<*>(k,k) 


k-1  ^ 

rxx(k)  - E <|>(k  - 1 , j )rxx( k - J) 



1 - E 4>(k  - 1,  j)r  (j) 
j = l XX 


k = 2,  3,  ....  n/4,  where, 


4>(k, j ) = <t>(k  - 1,  j)  - 4>(k,k)  <p(k  - 1,  k - j). 


(2.4.6) 


we  may  possibly  be  able  to  obtain  some  insight  into  the  identification  of  the 
stochastic  model  and  its  order  that  will  best  fit  the  data,  [6].  The  sample 
autocorrelation  function  of  an  autoregressive  process  of  order  p tails  off, 
while  its  sample  partial  autocorrelation  function  has  a cutoff  after  lag  p. 
Conversely,  the  sample  autocorrelation  function  of  a moving  average  process 
of  order  g has  a cutoff  after  lag  q,  while  its  sample  partial  autocorrelation 
function  tails  off.  If  both  the  sample  autocorrelations  and  partial  autocor- 
relations tail  off,  a mixed  process  is  suggested. 

2.4,2  The  Fitting  Procedure 

The  process  of  fitting  any  one  of  the  three  stationary  stochastic 
models  under  consideration  usually  involves  two  stages 

i)  deciding  the  order  of  the  process; 
ii)  estimating  the  appropriate  set  of  parameters. 

The  criterion  for  selecting  the  order  of  the  process  that  will  give  the  best 
fit  is  the  residual  variances  of  different  orders  of  the  process  fitted  to  the 
data.  To  compute  the  residual  variances,  it  is  necessary  to  estimate  the 
parameters  for  the  different  order  processes.  The  residual  variances  are  then 
plotted  against  the  order;  the  minimum  residual  variance  will  correspond  to 
the  correct  order  for  the  process.  After  this  has  been  completed  for  each 
process  (the  autoregressive,  the  moving  average,  and  the  mixed  autoregressive- 
moving  average),  we  compare  the  minimum  residual  variances;  the  minimal  one 
will  correspond  to  the  process  (and  its  order)  that  will  best  describe  the 
data.  This  procedure  is  quite  well  suited  for  use  on  a high-speed  computer. 
Note,  that  when  fitting  a model  to  a given  set  of  observations,  we  keep  in 
mind  the  principle  of  parsimony,  [6].* 

* Chapters  1 and  9 
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As  pointed  out  by  R.  A.  Fisher,  [7],  for  tests  of  goodness-of-fit  to 
be  relevant,  it  is  necessary  to  make  efficient  use  of  the  data  in  the  fitting 
process.  Hence,  to  obtain  efficient  estimates  of  the  parameters,  we  shall 
use  maximum  likelihood  estimates,  or  approximate  maximum  likelihood  estimates 
for  the  parameters  that  constitute  the  different  models.  The  asymptotic 
properties  of  maximum  likelihood  estimates  are  usually  derived  for  independ- 
ent observations,  but  as  was  shown  by  Whittle,  [8],  they  may  be  extended  to 
cover  stationary  time  series. 

Suppose  there  exists  a non-stationary  series,  xt,  t = 1,  2,  ...,  n +1 , 
generated  by  a mixed  autoregressive-moving  average  process  of  order  (p,  q), 
whose  first  difference,  yt,  t = 1,  2,  ...,  n,  is  stationary.  We  desire  to 
fit  a stationary  mixed  process  of  order  (p,  q)  to  the  y's;  that  is 

yt  = alyt-l  + + apyt>P  + Zt  " SlZt-1  ' ' Sq2t-q  . (2.4.7) 

Without  loss  of  generality,  one  can  assume  that  when  d > 0,  = 0.  We  can 

express  (2.4.7)  as 

zt  = ‘Vt-l  - •••  - Vt-P  + S1Zt-l  + •••  * Vt-q  • (2'4-8) 

A recursive  technique  can  now  be  used  to  obtain  the  conditonal  sum  of  squares 
function.  By  conditional  sum  of  squares  function,  we  mean  that  the  sum  of 
squares  function  is  conditional  on  the  starting  values  assigned  to  the  para- 
meters (a' s and  B's)and  for  the  y's  and  z's  previous  to  time  t.  This  simple 
numerical  technique  will  recursively  build  up  the  log-likelihood  function. 

Now,  assuming  the  Zt  process  is  normal  with  mean  zero  and  variance  a^, 
the  joint  probability  density  of  the  Z's  is 

f(zr  z2,  ....  zn)  - cr-n  exp[-  S zJ/2o*]  . (2.4.9) 

Given  a particular  set  of  data,  yt,  t = 1,  2,  ...,  n,  the  conditional  log- 

likelihood  associated  with  the  parameter  values  (a^ ap,  B-j , ....  Bq, 

az),  conditional  on  the  choice  of  starting  values  for  the  y's  and  z's,  is 
given  by 
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fc(a.|,  . .., 


dp.  6-]  > •••.  6p»  o2)  = -nixi  az 


and  the  conditional  sum  of  squares  function 


S(oti , . . • > otp » Bi , . • • > 6g ) 

2oz 

(2.4.10) 


S(ct-j  , ■ . . i cXp , Bi  , • ••>  8q) 


= E zj(a,,  ...,  a , 8,,  ...,  8„  (starting  values). 
t=l  U P 1 H 


(2.4.11) 


Notice  that  the  conditional  likelihood  (2.4.10)  involves  the  data  only 
through  the  conditional  sum  of  squares  function  (2.4.11).  It  follows  that 
contours  of  (2.4.10)  for  any  fixed  value  of  az  in  the  space  (a^  , ....  a , 

Bp  ...,  8 > a2)  are  contours  of  (2.4.11),  that  these  maximum  likelihood 
estimates  are  the  same  as  the  least  squares  estimates,  and  that,  in  general, 
one  can,  on  the  normal  assumption,  study  the  behavior  of  the  conditional 
likelihood  by  studying  the  conditional  sum  of  squares  function. 

Thus,  we  will  obtain  least  square  estimates  (maximum  likelihood 
estimates)  by  minimizing  the  sum  of  squares  function.  Note  that  the  para- 
meter values  are  selected  to  recursively  calculate  the  sum  of  squares 
function  such  that  they  will  satisfy  the  stationarity  and/or  the  inverti- 
bility  conditions  of  the  stationary  stochastic  model.  The  residual 
variances  can  then  be  obtained  by  dividing  the  sum  of  squares  function  by 
the  appropriate  degrees  of  freedom.  The  fitting  procedure  will  be  described 
in  greater  detail  for  the  autoregressive,  the  moving  average,  and  the  mixed 
autoregressive-moving  average  process  in  Section  3. 

2.4.3  The  Backward  Filter  and  Diagnostic  Checks 

Having  selected  the  stationary  stochastic  model  and  its  order  that  best 
describes  the  data  and  having  estimated  its  parameters,  diagnostic  checks  are 
conducted  on  the  model  to  determine  its  adequacy.  If,  in  the  identification 
stage  of  the  analysis,  a suitable  value  for  d was  found  to  be  different  from 
zero,  the  model  has  been  fitted  to  the  stationary  (filtered)  data,  not  to  the 
observed  non-stationary  data.  Hence,  to  use  the  fitted  model  to  forecast 
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future  values  of  the  observed  non-stationary  time  series,  we  introduce  the 
concept  of  a "backward  filter."  The  backward  filter  is  essentially  a simple, 
but  very  important , technique  that  allows  the  analysis  of  non-stationary  time 
series  using  the  well-established  methods  of  stationary  time  series. 

For  a better  understanding  of  the  above  concept,  consider  the  following 
example:  The  difference  equation, 

h 3 ’ B)d  xt  • 

was  used  to  filter  the  observed  non-stationary  time  series,  xfc,  t 3 1,  2, 

...»  n + d.  To  the  stationary  series,  yt,  t 3 1,  2,  ...,  n,  a 2nd  order 
autoregressive  model  is  fitted;  that  is 

/V  A 

yt  » ct]yt_1  + a2yt_2  + z%  . (2.4.12) 


To  transform  the  model  (2.4.12),  fitted  to  the  stationary  data,  to  the 
non-stationary  data,  xt,  we  simply  replace  yt  in  the  model  with  (1  - 8) 
that  is 


(1  - B)d  xt  • ^(1  - B)d  xt_1  ♦ a2(l  - B)d  xt_2  + zt  . 


For  d 3 1 and  simplifying,  we  have 


(2.4.13) 


xt  = *hxt-l  + *2xt-2  + 4>3xt-3  + zt  (2.4.14) 

A A A /V 

where  =1  + ot^ , <i>2  * a2  - , and  <j>3  = -a2>  The  $'s  are  linear  combina- 

tions of  the  a's  and  will  depend  upon  the  degree  of  differencing  needed  to 
filter  the  observed  series.  For  example,  if  d * 2,  we  have 


xt  3 'Vt-l  + 4>2xt-2  + ^t-S  + *4xt-4  + 2t  ‘ (2-4.15) 

A A A A A /V 

where  <t>-|  * 2 + , <J>2  3 a2  - 2ct^  - 1 , $3  3 - 2a2,  and  <)>^  * a2  . 

We  can  now  conduct  diagnostic  checks  on  the  fitted  model  (2.4.14)  to 
determine  the  goodness-of-fit.  Using  the  model  (2.4.14),  the  behavior  of  the 
observed  non-stationary  series  can  be  simulated.  By  plotting  the  simulated 
series  against  the  observed  series,  one  can  obtain  a visual  conception  of  the 
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goodness-of-fit.  More  substantially,  one  can  calculate  the  residuals 
incurred  by  subtracting  the  modeled  series  from  the  observed  series;  that  is 

A /S 

zt  * xt  - xt  . (2.4.16) 

If  the  model  is  adequate,  the  residuals  should  be  the  sum  of  the  z.  process 

A U 

plus  a factor  of  l//n;  that  is  as  n increases,  the  residual  z^'s  should 
behave  approximately  like  the  white  noise  2 t 1 s . Thus,  the  study  of  the 
residuals  could  indicate  the  existence  of  model  inadequacy,  in  particular, 
the  analysis  of  the  sample  autocorrelation  function  of  the  residuals,  [4  ]. 

If  the  form  of  the  model  was  known  to  be  accurate  and  the  true  parameter 
values  were  known,  then,  by  a result  of  Anderson,  [9  ],  the  sample  autocor- 
relations rzz( k)  of  the  z^'s  would  be  uncorrelated  and  distributed,  more  or 
less,  normally  about  zero  with  variance  1/n.  Therefore,  we  could  run  a 
statistical  test  to  determine  if  the  deviations  of  these  autocorrelations 
from  their  theoretical  zero  values  are  significant. 

However,  in  practice,  we  do  not  know  the  correct  form  of  the  model  nor 
the  true  parameter  values.  The  residuals  obtained  from  equation  (2.4.16) 

A 

will  be  estimates  of  the  residuals,  zt's,  not  the  zt's.  Hence,  the  accept- 
ance of  the  hypothesis  that  the  sample  autocorrelations  of  the  residuals 
constitute  a purely  random  process  on  the  assumption  of  a standard  error  of 
l/v'n  can  be  very  dangerous,  [10].*  Further,  it  is  shown  in  [6]^  that,  by 
using  l//n  as  the  standard  error  for  rzz(k),  the  statistical  significance  of 
the  deviations  from  zero  of  the  sample  autocorrelations  at  low  lags  will  be 
underestimated,  while  at  moderate  or  high  lags,  1/^n  will  give  an  acceptable 
estimate. 

Instead  of  considering  the  sample  autocorrelations  of  the  residuals 
separately,  an  indication  is  often  needed  of  whether  or  not  the  first  K auto- 
correlations of  the  residuals,  taken  as  a whole,  indicate  inadequacy  of  the 
fitted  model,  [4  ].  (The  value  of  K can  be  taken  to  be  n/10).  Now,  assume 

A 

that  we  have  calculated  the  residual  z t ' s and  their  first  K autocorrelations, 
rZ2 (k),  k = 1,  2,  ...,  K,  resulting  when  a fitted  mixed  autoregressive-moving 
average  model  of  order  (p,q)  is  used  to  model  the  observed  series.  Then,  it 
can  be  shown,  [6],  that  if  the  fitted  model  is  adequate, 

K „ 

Q * n £ rMk)  (2.4.17) 

k*l  zz 

★ 
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is  approximately  chi-square  distributed  with  (K-p-q)  degrees  of  freedom. 

Here,  n is  the  number  of  observations  used  to  fit  the  stationary  model.  If 
the  model  is  inadequate,  the  Q value  will  be  inflated.  Therefore,  one  can 
make  an  approximate,  general,  or  "postmanteau"  test  of  the  hypothesis  of 
model  inadequacy  with  the  information  available  by  referring  a calculated 
value  of  Q to  a table  of  percentage  points  of  the  chi-square  distribution, 
[6]. 

Note  that  for  fitted  autoregressive  models  of  order  p and  for  fitted 
moving  average  models  of  order  q,  we  would  have  a chi-square  with  (K-p)  and 
(K-q)  degrees  of  freedom,  respectively. 

2.4.4  Forecasting  and  Updating 

After  having  obtained  a model  to  describe  the  observed  time  series  and 
having  confirmed  its  adequacy,  we  desire  to  use  it  to  forecast  future  values 
of  the  observed  series.  We  shall  now  illustrate  how  the  fitted  model  may  be 
used  to  obtain  minimum  mean  square  error  forecasts.  We  would  like  to  fore- 
cast a value  xt+Jl,  £ = 1,  2,  ...,  L steps  ahead,  when  we  are  presently  at 
time  t.  That  is,  the  forecast  is  said  to  be  made  at  origin  t for  a lead  time 
l.  Of  course,  the  shorter  the  lead  time  i,  the  more  accurate  one  can  expect 
the  forecast  to  be.  Also,  the  spacing.  A,  of  the  data  is  of  importance  in 
forecasting.  That  is  to  say,  if  the  data  is  recorded  daily,  then  it  would 
be  unrealistic  to  attempt  to  forecast  weeks  In  advance. 

To  derive  the  minimum  mean  square  error  forecasts  for  any  lead  time  4, 

first  consider  the  general  mixed  autoregressive-moving  average  model  fitted 

to  the  stationary  series  yt;  that  is, 

h ’ Vt-i  * Vm  * •••  ♦ Vt-p  * zt  * 6izt-i  - ••• 

Using  the  backward  filter  (1  - B)dxt  ■ yt,  we  have 

xt  * *lxt-l  + 4>2xt-2  * + ^p+dXt-p-d  + zt  ‘ Sl2t-1  *■*  8qZt-q  ’ 

(2.4.19) 

or  in  terms  of  the  backward  shift  operator  B 

(1  - t,B  - 0282  - - *p+dBP+d)xt  « (1  - BjB  - B2B2  - ...  - SqBq)zt  . 

(2.4.20) 


Vt-q 


(2.4.18) 
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We  may  write  an  observation  xt+^ 
a difference  equation, 


generated  by  the  process  (2.4.19)  either  as 


(tn  ' ‘Vtn-l  + ^2xt+?,-2  + •••  + ^p+dxt+)i-p-d 


+ zt+2,  ‘ 6lzt+2-l  ■*•••“  Bqzt+H-q  ’ 


(2.4.21) 


or  as  an  infinite  weighted  sum  of  current  and  previous  shocks  z,, 

J 


't+fc 


00 

= Z I p.z 


j-0  j tn-j  ’ 


(2.4.22) 


where  ipQ  = 1 and  the  weights  may  be  obtained  by  equating  coefficients  in 


(1  - d>-,B  - 


yd8P+d)0  + ^8  + <P2B2+  ...)  = 
0 - S-,B  - ...  - B Bq); 


(2.4.23) 


or  as  an  infinite  weighted  sum  of  previous  observations,  plus  a random  shock, 


t+l 


= Z 

M 


1TjXt+A-j 


+ z 


t+i 


(2.4.24) 


Now,  using  equation  (2.4.22)  and  the  assumptions  of  the  model,  as 
discussed  previously,  it  will  be  shown  that  the  minimum  mean  square  error 
forecast  at  origin  t,  for  a lead  time  l,  is  the  conditional  expectation  of 

Xt+H  at  time  t;  that  1S’ 

xt(«.)  * EtCxt+£3  . (2.4.25) 


Suppose  at  time  t we  are  to  forecast  xt+Jl  with  a linear  function  of  current 
and  previous  observations  xt,  x^,  xt_2’  •••  • Then,  as  shown  above,  it 
will  also  be  a linear  function  of  current  and  previous  z^,  zt  zt_2»  ... 
Suppose  the  best  (minimum  mean  square  error)  forecast  is  given  by 


XtU)  - **lt  ♦ **  * 


= Z 


* 


\b 

jZt+£-j 


(2.4.26) 
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Now , from  equation  (2.4.22), 


lt+4 


x>  = jo  Vt«-j  • j,  'Mi 


4-1 


(2.4.27) 


' jfo  ^ - Vzt«-j  ■ 

Squaring  the  above  expression  and  taking  the  expected  value,  we  have 

4-1 


E[(x 


f»-  x‘(t)  )!]  = E[(jfo  Vtu-j  * ih  ‘*j  ' Vzf*-j)!] 


Since  E(z^Zj)  3 0,  for  i / j and  a*  for  i = j,  we  have 

E^xt+4  " = 2 ^jaz  + 2 - <Pj)2  al  (2.4.28) 

3=0  J J J z 

The  above  expression  is  minimized  by  setting^*3  $ . , j = z,  i + i , . . . 5 ®; 
this  implies: 


xt(i),Vt  + Wt.i  + - • 


(2.4.29) 


A 

Now,  it  is  required  to  show  that  x t ( £)  as  given  in  equation  (2.4.29)  is, 
in  fact,  the  conditional  expectation  at  time  t of  xf  Since 


Co  , for  j > t 
* 1 

J (_zj  ’ for  J 1 4 

from  equation  (2 .4. 22),  we  have 

00 

Et(xtu)  S Et^0  ^jZt+4-j^ 

= j»0  Et[Zt+4-j] 


jaE/jztn-j  • 


(2.4.30) 


(2.4.31) 
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A 


This  result  also  has  been  shown  in  [ 6 ]* 

For  example,  if  a mixed  autoregressive-moving  average  model  of  order 
(p,q)  was  fitted  to  the  dth  difference  series  of  the  observed  non-stationary 
series  x^.. 


xt  = Vt-l  + •••  + VdXt-p-d  + zt  ' Slzt-1  ' •••  - 6qzt-q  * 

(2.4.32) 

Now,  to  forecast  ahead  by  a lead  time  l,  we  replace  t with  t + £in  (2.4.32) 

xt+£  = <f>lxt+)l-l  + •••  + 4>p+dxt+i-p-d  + zt+£  " Slxt+£-1  " •**  ' 8azt+£-q  • 

(2.4.33) 

The  minimum  mean  square  error  forecast  will  be  given  by 

Et  ^Xt+J^  = ^1  Et^xt+«.-l + *•'  + ^p+dEt^xt+£-p-d^  + Et^zt+£^ 

"6lEt^zt+£-1-l  ”•  8qEt^zt+£-q-^  • ^ ^ 

As  discussed  previously,  the  conditional  expectation  can  be  obtained  from 


Et^xt+£^  = xt^  and  Et^zt+t^  " ® (2.4.35) 

for  i = 1 , 2,  . . . L,  and 

Et^xt->^  = xt-£  and  Et^zt-l^  3 zt-£  » (2.4.36) 

for  l - 0,  1,  2,  ...,  L.  Thus,  we  can  rewrite  equation  (2.4.35)  by  using  the 
following  rules: 

i)  The  data  points  x^+Jl  (£  * 1,  2,  ...»  L),  which  have  not  yet  been 
realized,  are  replaced  by  their  forecasts  *tU)  at  origin  t. 
ii)  The  errors  zt+?_  (i  = 1,  2,  ....  L),  which  cannot  be  calculated  until 
xt+£  1S  rea^1zed’  are  replaced  by  their  unconditional  expectation 
of  zero.  (See  section  3.3.4). 

We  can  also  rewrite  equation  (2.4.36)  by  the  following  rules: 

i)  The  data  points  xt_^  [l  * 0,  1,  2,  ....  L),  which  have  already  been 
realized  at  orig-'  t,  assume  their  realized  value, 
ii)  The  errors  zt_^  (?.  = 0,  1,  2,  ...  L),  which  have  happened  at  origin 
t,  are  calculated  from  xt_^  - xt,£_-|(l). 


* Chapter  5 


24 


The  variance  of  the  i steps  ahead  forecast  error  for  any  origin  t is  the 
expected  value  of 


e|U)  = [xt+r  xtU)]2  . (2.4.37) 


It  has  been  shown,  [12],  that  the  variance  of  the  lead  time  i is  given  by 


£-1 

var(£)  = [1  + L 0?]  a2  (2.4.38) 

j-1  J 


where  the  weights  0.  are  given  by 
J 


9j  = 0, 
9o  = 1 


j < 0 


0^  = <j>1  - ^ 

A 

^2  = ^1®1  + ^2  ” 92 


(2.4.39) 


9j  = *l0j-l  + •••  + Vd9j-p-d  - Bj’  J = 1.  2, 


For  j greater  than  q and  p+d-1*  equation  (2.4.39)  can  be  reduced  to: 


9j 


Vm 


+ »28j-2  + 


"p+d^j-p-d 


(2.4.40) 


Note  that  when  one  has  an  autoregressive  model,  the  B's  are  zero;  similarly, 
when  one  has  a moving  average  model,  the  <t>'s  are  zero. 

We  may  express  the  accuracy  of  the  forecasts  by  calculating  probability 
limits  on  each  forecast.  The  probability  limits  are  such  that  when  the 
realized  value  of  the  time  series  occurs,  it  will  be  included  within  these 
limits  with  the  stated  probability.  An  estimate  s2  of  the  variance  a2  in 
(2.4.38)  is  obtained  from  the  time  series  data,  s2  will  be  the  residual  sum 
of  squares  obtained  in  the  fitting  procedure  divided  by  the  number  of  obser- 
vations used  in  calculating  it.  When  the  number  of  observations  is,  say,  at 
least  50,  one  can  approximate  the  (1-a)  probability  limits  as, 

* That  is,  j > p+d-1  if  p+d-1  > q;  j > q if  p+d-1  < q 
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(2.4.41) 


xt+t  = xt^  ± taf2  ^ + ^ 0P  1 sz  » 

where  t^  is  the  deviate  from  the  student's  t-distribution. 

We  are  usually  interested  in  forecasting  values  of  an  observed  series 
for  z = 1 , 2,  ....  L lead  times  in  the  future  from  some  origin  t.  Once  these 
forecasts  are  obtained  and  a new  piece  of  information  is  realized,  we  may 
adjust  or  update  the  original  forecasted  values.  The  new  forecast  will  be 
related  to  the  old  by 


Xt+1  U)  = + 1 )+9!lZt+l 


(2.4.42) 


From  the  above  equation,  we  see  that  the  t-origin  forecast  of  xt+j^  may  be 
updated  to  become  the  t + 1 origin  forecast  of  the  same  xt+^+-j , by  adding  a 
constant  multiple  of  the  one-step-ahead  forecast  error  zt+^ , where  zt+^  is 
given  by 

/v 

Zt+1  = xt+i  - xt(l)  , (2.4.43) 

and  the  multiplier  9^  is  given  by  (2.4.39)  and  (2.4.40) 


2.5  THE  SPECTRUM 


After  fitting  the  appropriate  model  to  the  data,  additional  information 
can  be  obtained  from  the  filtered  (stationary)  series  concerning  the  distri- 
bution of  the  variance  with  respect  to  frequency.  The  Fourier  transform  of 
the  autocovariance  function  is  another  means  by  which  a stationary  stochastic 
process  can  be  characterized.  One  is  forced  to  conclude,  however,  that 
Fourier  analysis  breaks  down  when  it  is  applied  to  non-stationary  time  series. 
The  reasons  are  rather  obvious;  the  theory  behind  Fourier  analysis  is  based 
on  the  assumptions  that  the  amplitudes  are  fixed,  as  well  as  the  frequencies 
and  phases.  This  is  not  the  case  when  we  have  to  deal  with  time  series. 

Random  changes  in  frequencies,  amplitudes,  and  phases  are  found  here  due  to 
the  nature  of  time  series. 

The  sample  spectrum  is  the  Fourier  transform  of  the  sample  autocovari- 
ance function.  Its  graph  shows  how  the  variance  of  the  realization  of  a 
stochastic  process  is  distributed  with  respect  to  frequency.  Furthermore, 
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having  the  sample  spectrum,  one  can  obtain  the  sample  autocovariance  function 
by  taking  the  inverse  Fourier  transform  of  the  sample  spectrum.  That  is,  the 
sample  spectrum  of  cxx(f)  is  given  by 


Cxx(f)  = /-Tcxx(u)e'j27rfUdu’  • °°<f 


< 00 


and  the  sample  autocovariance  function  of  Cx  (f)  is  given  by 


(2.5.1) 


cxx(u)  = Oxx(0^df>  -T  < u < T 


(2.5.2) 


for  a continuous  time  series  x(t),  while  for  the  discrete  case, 

Cxx(f)  - A V _.c,Y(k)e'j27rfkA  , - <f<l 


k=- ( N- 1 ) 


xx 


2a  - J ^IK 


(2.5.3) 


and 


1 


j2Trfu 


cxx(u)  = f IK  CYY(f )eJiTTTUdf , -NA  < u < NA 


1_ 

2A 


xx 


(2.5.4) 


The  theoretical  spectrum  is  defined  by  taking  the  limit,  as  the  period 
T tends  to  infinity,  of  the  expected  value  of  the  sample  spectrum.  That  is. 


r (f ) = 11m  E[C  (f)] 
T-x» 


■ r.  Yxx(u)e-J2”fuciu 


(2.5.5) 


Similarly,  the  theoretical  autocovariance  function  of  the  continuous  time 

\ 

series  x(t)  is  obtained  by  taking  the  inverse  Fourier  transform  of  the  theore- 
tical spectrum.  That  is, 


(2.5.6) 


The  graph  of  the  spectrum  rxx(f),  as  a function  of  frequency,  shows  how 
the  variance  of  the  x(t)  process  is  distributed  with  respect  to  frequency. 
Similarly,  the  graph  of  the  sample  spectrum  C (f)  shows  how  the  variance  of 
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a realization  of  a stochastic  process  is  distributed  with  respect  to 
frequency  over  a period  of  length  T. 

Quite  frequently  in  practice  we  have  to  compare  time  series  with  differ- 
ent scales  of  measurements.  Thus,  it  is  necessary  to  normalize  the  spectrum. 
Normalization  in  this  case  is  done  simply  by  dividing  the  theoretical  spec- 
trum given  by  equation  (2.5.5)  by  the  variance  of  the  process  a*.  That  is, 

rxx(f) 

Wf>  ■ -*%-•  (2-5-7) 

where 

= Yxx<°>  * '-rxx(f)df  ■ 

The  expression  given  by  equation  (2.5.7),  £ (f),  is  called  the  spectral 

density  function.  The  function  £ ( f)  represents  the  Fourier  transform  of 

)\A 

the  autocorrelation  function  due  to  the  relationship  between  the  autocovari- 
ance and  the  autocorrelation  function.  The  spectral  density  function  is 
non-negative  and  integrates  to  unity,  thus  resembling  the  definition  of  a 
probability  density  function. 

2.5.1  Properties  of  the  Spectrum 

Sometimes  it  is  highly  desirable  to  obtain  the  spectrum  of  the  output 
from  a linear  system  when  the  input  is  a stationary  process.  In  the  present 
study,  we  will  be  primarily  interested  in  the  case  where  the  input  is  white 
noise.  The  general  rule  is  given  in  Jenkins  and  Watts  [3]*,  which  states  that 
"the  spectrum  of  the  output  from  a linear  system  is  obtained  from  the  spectrum 
of  the  input  by  multiplying  by  the  square  of  the  frequency  response  function." 
Thus,  in  our  case,  the  input  to  the  system  will  be  white  noise  denoted  by 
z(t)  whose  spectrum  is  given  by  T (f)  = a*,  and  the  output  process  x(t)  is 
a linear  process  and  its  spectrum  is  given  by 

rxx(f)  = CT|lH(f)  I2 , -00  < f < 00  (2.5.8) 

* Chapter  6 
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or 


rxx(f)  = lH(f)i2rzz(f)»  -k-f-k  (2.5.9) 


where  H(f)  is  the  frequency  response  function. 

In  our  presentation  of  the  spectral  analysis  of  the  ionospheric, 
man/machine  interface, and  climatological  information,  we  shall  utilize  the 
spectrum  of  the  three  basic  models  we  discussed  in  sections  2.1  through  2.4, 
namely,  the  autoregressive,  the  moving  averages,  and  a mixture  of  auto- 
regressive-moving averages.  Thus,  in  what  follows,  we  shall  give  the  basic 
definitions  of  the  above  models  for  both  the  continuous  and  the  discrete 
case. 

i )  Continuous  First  Order  Autoregressive  Process 
The  spectrum  of  a continuous  realization  of  a stochastic  process  x(t)  is 
given  by 


rxx(f)  = , 

xx  1+  ( 2irf  T ) 2 


f < 


(2.5.10) 


Note  that  from  the  form  of  the  theoretical  spectrum,  r (f),  one  can  conclude 
that  most  of  the  power  (variance)  is  concentrated  at  low  frequencies. 

i i )  Discrete  First  Order  Autoregressive  Process 
The  spectrum  of  a discrete  realization  of  a stochastic  process  x(t)  is 
defined  by 


rxx(f> 


4C2 


l+a|-2a-|Cos  2rfA 


1 

IK 


< f < 


1 

IK  • 


(2.5.11) 


It  should  be  mentioned  that  if  is  negative,  the  spectrum  has  most  of  its 
power  concentrated  at  higher  frequencies,  while  when  is  positive,  the 
power  will  be  concentrated  at  lower  frequencies. 

i i i )  Continuous  Second  Order  Autoregressive  Process 
The  spectrum  of  a continuous  realization  of  a stochastic  process  x(t)  is 
given  by 


rxx<f> 


(a0-a24ir2f2): 


(2irfa1 ) 


< f < 


(2.5.12) 
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Here  the  distribution  of  the  power  of  the  process  will  depend  on  how  aQ,  a] 
and  a2  are  inter-related.  This  relationship  is  inherent  in  the  form  of  the 
characteristic  equation  of  the  process. 

i v )  Discrete  Second  Order  Autoregressive  Process 
The  spectrum  of  the  continuous  realization  of  a stochastic  process  x(t)  is 
given  by 


A at 


rxx<f> 


l+a^+a^"^-]  (l"a2)cos  2irfA-2a2cos  4irfA 


ie  i f ik 


(2.5.13) 

Again,  the  location  of  where  most  of  the  power  is  concentrated  is  character- 
ized by  the  values  of  and  a2. 

v)  General  Autoregressive-Moving  Average  of  a Continuous  Process 
The  spectrum  of  a continuous  realization  of  a stochastic  process  x(t)  is 
given  by 


r 


xx 


i 


bQ+blj27Tf+...+b!l(j2TTf)& 
a0+alj2*f+-  • •+a|T)(j27rf  )m 


-00  < f < oo  , 


(2.5.14) 


vi )  General  Autoregressive-Moving  Average  of  a Discrete  Process 
The  spectrum  of  a continuous  realization  of  a stochastic  process  (x(t)  is 
defined  by 


rxx<f> 


l+S1e'J'2TrfA+. . .+g^e"';'27TfA£ 

1-cue-j27TfA-...-a  e-j2TTfArn 
1 jn 


< f < 


1 

IE 


(2.5.15) 


In  the  last  two  definitions,  various  peaks  or  spikes  will  appear  in  the 
spectrum  if  the  roots  of  the  corresponding  characteristic  equations  are 
complex.  The  behavior  of  the  above  will  be  displayed  in  Section  4. 


2.6  ESTIMATE  OF  THE  SPECTRUM 

In  the  previous  section,  we  have  given  a brief  discussion  of  the  basic 
definitions  and  properties  of  the  theoretical  spectrum.  With  respect  to  the 
aims  of  the  present  study,  we  shall  give  in  this  section  a smoothed  estimate 
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of  the  theoretical  spectrum  and  some  of  its  properties  that  are  essential  in 
the  final  analysis  of  the  experimental  data. 

Bartlett,  [ll],  introduced  a very  useful  procedure  for  estimating  the 
spectrum  of  a stochastic  realization.  This  estimate  of  the  spectrum  for  a 
purely  random  process  is  given  by 


. n-1  n-1 

C„(f)  = [(  Z ztcos  ZirftA}2  + { E z^sin  2irfA}2] 

t=-n  t=-n 


= j|  [A2(f ) + B2(f)]  , - £ f <.  , 

where 

n-1 

A2(f ) = { E ztcos  2irftA}2 
t=-n 

and 

n-1 

B2(f ) = { I z.sin  27ifA}2  . (2.6.3) 

t=-n 


(2.6.1) 


(2.6.2) 


Bartlett's  procedure  consists  of  splitting  up  the  series  of  length  N 
into  k sub-series  of  length  p evaluates  a sample  spectrum  C ^(f)  for  each 
of  the  k sub-series,  i = 1,  2,  ...,  k,  and  finally  takes  the  mean  of  the  sub- 
series as  his  estimator  at  frequency  f.  That  is, 

Zzz(f)  = \ Z czz{1)(f)  • (2-6.4) 

The  estimate  given  by  equation  (2.6.4)  is  called  a smoothed  spectral  estimate 
at  frequency  f and  the  method  to  obtain  (2.6.4),  Bartlett's  smoothing 
procedure.  More  generally,  Bartlett's  smoothing  procedure  suggests  that  a 
smoothed  spectral  estimator  of  a stochastic  realization  x(t)  is  given  by 

Cxx(f)  * 0(u)cxx(u)^du  = rjxx(u)e‘j27rfudu  . (2.6.5) 

The  smoothed  sample  spectrum  will  have  a smaller  variance  than  the 
unsmoothed  sample  spectral  estimator,  Cxx(f).  The  smoothed  estimate  of  the 
spectrum  is  a function  of  the  type  of  lag  window,  w(u),  one  utilizes. 
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However,  in  all  cases  the  following  properties  hold  with  respect  to  w(u): 

i)  w(0)  = 1; 

ii)  w(u)  = w(-u); 

iii ) w(u)  =0,  | u | > T; 

iv)  w(u)=0,  |u(  > M,  M<T. 

In  section  2.8  the  above  concept  of  the  window  will  be  explored 
extensively.  As  will  be  stated,  the  bandwidth  of  a spectral  window  will  be 
defined  as  follows: 


b = 


1 

rymdf 


(2.6.6) 


Another  form  of  the  bandwidth  is  the  standardized  bandwidth,  b^ , which 
is  given  by  placing 


and 


b 


02(f)df 


We  can  conclude  that  the  variance  of  a spectral  estimator  is  inversely 
proportional  to  the  bandwidth  of  the  spectral  window.  Also  the  degrees  of 
freedom,  v,  of  the  smoothed  estimator  are  directly  proportional  to  the  band- 
width of  the  spectral  window  due  to  the  relationship 


v 


2T 

/>2(u)du 


2(J)b, 


(2.6.7) 


On  the  other  hand,  the  bias  is  directly  proportional  to  the  bandwidth  of  the 
window. 

The  notion  of  the  smoothed,  spectral  density  estimate  denoted  by  IT  (f) 
shall  now  be  introduced.  This  estimate  is  defined  as  follows: 
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l-l 

(f)  = 2 {1  + 2 Z r (k)w(k)cos  2irfk}, 
k=l  xx 


0<  f < 


1 

2 ’ 


where, 


00 

m 


The  mean  smoothed  spectral  density  r (f) 

~T2 


is  given  by 


(2.6.8) 


(2.6.9) 


rx*(f>  L'1 

= 2 [1  + Z p (k)w(k)cos  2irfk]  . (2.6.10) 

a x k=l  xx 

Equation  (2.6.10)  is  the  expected  value  of  the  smoothed  spectral  denisty 
estimator.  By  plotting  rxx(f)  versus  "r  (f)  , as  shall  be  done  with  the 


ionospheric  phenomenon  under  investigation,  one  will  be  able  to  detect  how 
the  bias  varies  with  frequency. 

Also,  the  variation  of  bias  with  bandwidth  can  be  observed  by  simply 
plotting  different  curves  for  the  different  values  of  L,  the  truncation 
point.  One  other  relationship  can  be  plotted  as  well.  This  plot  consists 
of  R (f)  versus  T (f)  and  displays  how  the  variance  varies  with  respect  to 

AA  AA 

frequency.  a x 

Similarly,  the  variation  of  variance  wi.th  bandwidth  can  be  observed 
simply  by  plotting  the  estimate  with  different  values  of  L. 

Jenkins  and  Watts,  in  their  efforts  to  point  out  the  criteria  for 
determining  optimal  lag  windows,  argue  that  window  carpentry  is  not  as 
important  as  window  closing.  In  Section  4 of  the  present  study,  we  will 
examine  these  factors  in  detail  with  actual  data  to  determine  an  optimal  lag 
window. 

Two  criteria,  which  would  seem  logical  at  first  sight  but  whose  value 
is  debatable  in  determining  an  optimal  lag  window,  are  listed  as  follows. 
(They  are  classified  as  the  "optimality  approach  to  smoothing."): 
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i)  The  mean  squared  error  criterion, 

E^*x<f>  - rxx<f))!^  • 
and 

ii)  The  integrated  mean  squared  error, 

C.  Bt^xx<f>  - I'xx(f)>!ldf  ■ 

Nevertheless,  the  above  optimality  criteria  are  rejected,  a fact  which 
brings  into  play  the  formulation  of  high  stability  and  high  fidelity  as 
useful  criteria.  They  are  widely  accepted  in  practice.  The  previous 
criteria,  namely,  the  mean  squared  error  and  the  integrated  squared  error, 
are  not  useful  because  they  are  arbitrary;  they  have  no  flexibility  (strictly 
mathematical);  they  do  not  allow  an  a priori  design  and  analysis  of  the  data; 
and  they  only  indicate  what  is  best  on  the  average.  Therefore,  the  two  main 
requirements  for  estimating  the  theoretical  spectrum  r (f)  as  accurately  as 
possible  are: 

i)  High  fidelity,  which  implies  that 

rxx<f>  - rxx<f>  * B<f> 

be  small,  and 

ii)  High  stability,  which  implies  that 
Var^f)]  : ({1  ) 

be  small . 

High  fidelity  and  high  stability  are  two  conflicting  requirements.  In 
minimizing  the  covariance,  we  increase  the  bias;  and,  conversely,  by  mini- 
mizing the  bias,  we  increase  the  variance.  An  ideal  situation  would  be  one 
where  M is  large  enough  for  high  fidelity  and  M is  small  enough  for  high 
stability.  The  logic  here  dictates  a compromise  of  some  form. 

In  summary,  we  can  say  that  smoothing  the  estimate  of  the  theoretical 
spectrum  consists  of  determining  the  shape  or  the  mathematical  form  of  the 
window  (section  2.8  discusses  window  carpentry)  on  one  hand  and  the  value  of 
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the  bandwidth  (window  closing)  on  the  other.  By  window  closing  we  mean 
computing  window  estimates  using  a wide  bandwidth  and  then  progressively 
using  narrower  bandwidths  until  we  achieve  a state  of  high  fidelity  and  high 
stability.  The  numerous  and  various  problems  involved  will  be  studied  in 
subsequent  chapters. 


2.7  THE  CROSS  SPECTRUM 


In  this  section,  the  concepts  we  have  already  dealt  with  shall  be 
extended  so  as  to  be  able  to  treat  pairs  of  time  series  x-j(t)  and  x2(t). 
The  autocovariance  function  of  the  stochastic  realization  is  given  by 

Yx  x (“)  * E[(X1(t)-U1 ) (X1(t+u)-u1)]  = Y-| i , (2.7.1) 

where  the  expression  for  can  be  obtained  by  changing  the  subscripts, 
and  the  cross-covariance  function  is  given  by: 

Yxx  (u)  » EC(X1(t)-y1)(X2(t+u)-u2)]  = Y 1 2 (2-7.2) 


where 


ui  = E[X.j (t)] , i = 1,  2. 

The  cross  correlation  function  is  given  by 

Y12(u)  Y12(u) 

Pl2<U>  ' /;'ToT”ToT  ' ‘W 

/ Y-j  i Y22 

The  cross  covariance  function  is  estimated  by 


(2.7.3) 


1 

T 


T/2-u 

(X,(t)-T,)  (X-(t+u)-I9)dt,  0 < u < T 
-T/2  1 1 c C 


xlx2 


(u) 


1 

T 


T/2 


-T/2+u 


(X] ( t)-X1 ) (X2(t+u)-X2)dt,  -T  < u < 0 


(2.7.4) 


35 


where 


X^tjdt, 


1-1,2 


(2.7.5) 


An  estimate  of  the  cross  correlation  function  is  obtained  by 


r„  „ (k)  - 


c w (k) 


xlxl  x2x2 


(2.7.6) 


Using  the  above  equations  as  a foundation,  one  can  define  the  theoretical 
autospectrum,  T-j^f),  and  sample  autospectrum,  C^(f),  by 


rT|(f)  = Oll(u)e*jMUdu  ’ 


(2.7.7) 


C-j  i ( f ) = /Tc11(u)e'j2lTfudu, 


(2.7.8) 


respecti vely.  Hence,  we  arrive  at  the  notion  of  the  cross  spectrum.  As  in 
the  univariate  case,  the  sample  cross  spectrum  is  the  Fourier  transform  of 
the  sample  cross  covariance  function.  That  is. 


C12(f)  * /ITc12(u)e"j2Trfudu 


(2.7.9) 


Note  that  the  inverse  Fourier  transform  of  C^ff)  gives  rise  to  c-^f) 
Another  form  of  C^(f)  that  i?  commonly  used  is 


C12(f)  - A12(f)e 


JFl2(t)  ^ l 


12(f)  - JQ12(f) 


(2.7.10) 


which  is  the  product  of  a real  function  A^(f),  called  the  sample  cross 
amplitude  spectrum , and  a complex  function  called  the  sample  phase  svectrum, 
jF12(f)  . 
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The  theoretical  cross  spectrum  is  the  Fourier  transform  of  the 
theoretical  cross-covariance  function;  that  is, 

ri2(f)  = Ol2(u)e‘j2irfUdu  = 

= A12(f)  * j4,12(f)  • (2.7.11) 

which  is  the  product  of  the  cross  amplitude  spectrum  and  the  phase  spectrum. 
The  expected  value  of  the  cross  spectrum  estimate  is  given  by 


E[CX  X (f)]  = ^T(1  - 4J'hx1x2(u)e"j27rfUdu  • (2-7.12) 

As  T goes  to  infinity,  the  above  mean  value  approaches  the  cross  spectrum. 
Therefore, 

lim  E[C  (f)]  = r (f)  = /”  yx  x (u)e'j27rfudu  < f < • . 

T-o,  xlx2  xlx2  -00  xix2 

(2.7.13) 

Furthermore,  the  co-spectrum , A-^f),  is  defined  as 
A12(f)  = 012(u)e'j2irfUd“ 


^ /^0{y12(u)  + Y-J 2^ “u ^ c°s)  2irfudu  , 


where 


*12(u)  3 2^Yi2(u)  Y-|2^“u ^ * 

Similarly,  the  sample  co-spectrum , L^f)*  i$  defined  by 

L12(f)  - /V12(u)e’j2lTfudu 


(2.7.14) 


(2.7.15) 


3 2 T_y{c-|2(u)  - C^2(-u)}COS  2irfudu 


(2.7.16) 


'•V 


where 


^12^)  - 2^ci2^u^  ” c12^"u) ^ " 

On  the  other  hand,  the  quadrature  spectrum,  ^i^f),  is  defined 

V12(f)  = 012(u)e"j2lTfudu 
1 00 

= 2" /_„{ Y12(u)  - Y12(-u)  sin  }2irfudu  , 


where 


tp12(u)  = \ ^12^)  ~ Y12(-u)}  » 
and  the  sample  quadrature  spectrum,  Q.j2(f),  is  defined  by 

Q12(f)  = /Tq12(u)e‘j27rfudu 


where 


= j/^j{c12(u)  - c12(-u)}sin  Trfudu 


^12^U)  ""  2^c-| 2 ( u ) " c]2^’u^ 


Thus,  with  the  above  definitions,  we  can  define  the  cross 
spectrum  as  follows: 


c12(f)  * [r12(f)|  * /l^Trrnf^TFT  . 

Also,  the  sample  cross  amplitude  spectrum  is  defined  by 

A12(f)  * |C12(f)l  = / A|2(fr+  • 


(2.7.17) 
by 

t 

i 

(2.7.18) 

(2.7.19) 

(2.7.20) 

amplitude 

(2.7.21) 

(2.7.22) 
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On  the  other  hand,  the  phase  spectrum  can  be  written  by 


$12(f)  = arctan  , (2.7.23) 

and  the  sample  phase  spectrum  is  given  by 

-Q12(f) 

F12(f)  = arctan  • (2.7.24) 

Having  a graph  of  the  cross  amplitude  spectrum,  one  can  detect  whether 
frequency  components  in  one  series  are  associated  with  large  or  small  ampli- 
tudes at  the  same  frequency  in  the  other  series.  The  graph  of  the  phase 
spectrum  helps  to  determine  whether  frequency  components  in  one  series  are 
in  phase  or  out  of  phase  (lag  or  lead)  with  the  components  at  the  same 
frequency  in  the  other  series. 

The  cross  amplitude  spectrum  and  the  phase  spectrum  would  suffice  to 
provide  a complete  description  of  a bivariate  stochastic  process.  However, 
a more  efficient  spectrum,  namely,  the  coherency  spectrum,  will  be  introduced 
in  sub-section  2.7.1  to  take  the  place  of  the  cross  amplitude  spectrum. 

For  the  discrete  case,  we  simply  replace  the  integral  with  a sum  and 
make  the  necessary  notational  changes. 


The  squared  coherency  spectrum  is  the  plot  of  k|2(f)  versus  frequency. 
The  cross  amplitude  spectrum  o-j2(f)  is  a measure  of  the  covariance  between 
the  two  time  series  ( t)  and  x2(t)  at  frequency  f.  F^(f)  is  the  variance 
of  the  input  at  frequency  f,  and  G(f)  is  the  gain  of  the  spectrum  defined  by 
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G(f)  » 


(2.7.26) 


a12(f) 

= r7T^  • 


In  general,  the  coherency  spectrum  plays  the  role  of  a correlation 
coefficient  with  respect  to  the  frequencies.  When  r (f)  = 0,  the  squared 
coherency  is  equal  to  one.  When  the  output  is  nothing  else  but  noise,  the 
squared  coherency  is  equal  to  zero. 

The  usefulness  of  the  squared  coherency  spectrum  lies  in  the  fact  that 
dimensions  do  not  enter  the  picture  when  the  correlation  is  measured  with 
respect  to  frequency.  Unlike  the  squared  coherency  spectrum,  the  cross 
amplitude  spectrum  depends  on  the  dimensions  of  x-j(t)  and  x2(t).  This  is  the 
reason  why  the  squared  coherency  spectrum  is  preferred  over  the  cross  ampli- 
tude spectrum,  and,  together  with  the  phase  spectrum,  it  gives  us  a complete 
picture  of  the  cross  correlation  properties  of  two  time  series. 


2.8  THE  ROLE  OF  THE  LAG  WINDOWS 


One  of  the  goals  of  the  present  study  is  to  perform  spectral  analysis 
on  univariate  and  bivariate  stochastic  realizations  obtained  from  ionospheric 
soundings.  The  purpose  of  such  an  analysis  is  important  to  the  systems 
design  engineer  and  to  the  communication/ADP  scientist.  One  of  the  basic  and 
essential  factors  which  enters  in  the  mathematical  formulation  of  the  power 
spectrum  is  the  lag  window.  There  are  specifically  four  lag  windows  that  are 
commonly  used  in  spectral  analysis.  These  lag  windows  are  as  follows: 

i)  rectangular  window 

ii)  Bartlett's  lag  window 

iii)  Tukey's  lag  window 

iv)  Parzen's  lag  window. 


t 

( 


The  object  of  this  section  is  to  briefly  introduce  these  lag  windows. 


A specific  discussion  about  their  behavior  in  estimating  the  spectrum  will 


be  given  in  Section  4,  where  we  consider  the  analysis  of  vertical  incidence 
and  oblique  incidence  ionospheric  information. 

An  exact  determination  of  the  autocovariance  function  or  the  power 
spectrum  function  is  quite  impractical  since  it  would  require  both  a collec- 
tion of  pieces  of  infinite  length  and  an  infinite  number  of  computations. 

I 
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An  approximate  determination  is  therefore  proposed  in  order  for  workable 
estimates  to  be  obtained.  The  results  obtained  by  Tukey,  [12],  show  that, 
although  the  autocovariance  function  and  the  power  spectrum  are  Fourier 
transforms  of  each  other,  the  latter,  in  most  practical  situations,  is 
preferred  as  yielding  better  results.  By  reducing  the  data,  the  estimates 
will  be  subject  to  the  usual  sampling  variations  and  statistical  biases.  It 
should  be  pointed  out  again  that  estimates  of  the  power  spectrum  exhibit  bias 
and  sampling  variability  characteristics  much  easier  to  study  than  estimates 
of  the  autocovariance  function. 

The  unsmoothed  sample  spectrum  is  given  by 

Cxx(f)  - /Tcxx(u)e*j2lrfudu,,-»  < f <-  . (2.8.1) 

The  first  moment  of  C (f),  the  unsmoothed  sample  spectrum  estimator,  is 
given  by 


E[Cxx(f)]  = /ITE[cxx(u)]e'j27rfudu  , (2.8.2) 


which  can  be  written  as 

E[Cxx(f)]  = /ITYxx(u)(l  - . (2.8.3) 


using  the  fact  that 


E[cxx(u)] 


(2.8.4) 


Therefore,  using  the  convolution  theorem,  the  expected  value  of  the  unsmoothed 
sample  spectrum  can  be  written  as 

E[Cxx(f)]  = 0‘{S~-y  V rxx(f-g)dg  . (2.8.5) 

g 

In  expression  (2.8.5)  the  quantity 


T{- 


sin  ir  T 


*■} 


2 
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in  the  integrand  is  referred  to  as  the  spectral  window  of  the  sample  spectrum. 
The  Fourier  transform  of  the  spectral  window  is  the  lag  window  which  will  be 
denoted  by  w(u). 

In  actual  practice  we  utilize  a smoothed  form  of  expression  (2.8.5)  to 
obtain  an  estimate  of  the  theoretical  spectrum  given  by 

L-l  l 

Cxx(f)  * 2Ccxx(0)  + 2 Z cxx(u)w(u)cov  2tt  fk],  0 < f < j , 

(2.8.6) 

where  cxx(k)  is  the  sample  autocovariance  at  lag  k,  w(k)  being  the  lag 
window,  and  L the  truncation  point  of  the  series. 

It  is  clear  from  a practical  point  of  view  that  the  resulting  estimate 
of  the  power  spectrum,  equation  (2.3.6)  depends  on  the  choice  of  the  lag 
window  that  we  are  utilizing  in  our  estimate.  In  what  follows,  these  lag 
windows  shall  be  defined  and,  for  a more  thorough  investigation  of  their 
origin  and  properties,  the  recent  books  of  Jenkins  and  Watts,  [3],  and  Box 
and  Jenkins,  [ 6 1 » are  recommended. 

2.9  USEFUL  LAG  WINDOWS 


In  this  section  the  basic  useful  windows  mentioned  above  shall  be 
defined,  namely,  the  rectangular , Bartlett's , Tukey's , and  Parzen's  lag 
windows.  A specific  comparison  of  these  lag  windows  in  actual  problems  will 
be  given  in  Section  4. 

We  shall  denote  W(f)asthe  portion  of  the  integrand  given  by 
T{sl^f-Tf.}2  of  equation  (2.8.5).  That  is, 

W(f)  « T{s-lj?-y V . (2.9.1) 

Equation  (2.9.1)  is  called  the  spectral  window.  Its  Fourier  transform  is 
called  the  lag  window. 


2.9.1  The  Rectangular  Lag  Window 


Definition  2.9.1  The  function  given  by 


wR(u)  = 


1,  I u I < M 


0,  otherwise 


(2.9.2) 


is  called  the  rectangular-  lag  window.  The  Fourier  transform  of  wR(u)  is 


wR(f)  = - < f < - , 


(2.9.3) 


and  is  called  the  rectangular  spectral  window.  The  above  lag  window  was  the 
first  window  that  was  used  by  many  scientists,  especially  engineers.  However, 
due  to  its  mathematical  simplicity  in  being  able  to  characterize  complicated 
phenomena,  it  is  not  very  useful.  Further  details  to  this  addendum  will  be 
given  in  a later  discussion. 


2.9.2  Bartlett's  Lag  Window 


Definition  2.9.2  The  function  defined  by  the  expression 


wB(u)  = 


1 - ^ * u £ M 


(2.9.4) 


otherwise 


is  called  Bartlett's  lag  window.  The  Fourier  transform  of  Wg(u)  is  given  by 

«B(f)  - (2.9.5) 

and  is  called  Bartlett’s  spectral  \ window . The  above  lag  window,  which  was 
introduced  in  the  1950' s,  has  been  used  extensively  in  spectral  analysis. 

It  possesses  some  interesting  features  that  will  be  discussed  later.  However, 
its  side  lobes  are  much  larger  than  any  of  the  other  windows  known. 
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2.9.3  Tuke.y's  Lag  Window 


Definition  2.9.3  The  function  Wy,  defined  by  the  following  expression. 


Wy(u)  = 


jO  + cos  • M £ M 


(2.9.6) 


, otherwise. 


is  called  Tukey's  lag  window.  The  Fourier  transform  of  equation  (2.9.6)  is 
called  Tukey's  spectral  window  and  it  is  given  by 

wT(f>  - * i sin2TW(f;?M) 


2irM(f  + j M) 


u/sin  2irfM\  ,1  ^ 

+ M(  "2*fM"  } (T-T2fM]^  ’ -°°  1 f 1 00  • 


(2.9.7) 


Tukey's  window  possesses  the  property  of  having  most  of  its  power  concen- 
trated at  low  frequencies.  Furthermore,  Tukey's  window  results  in  smaller 
bias  in  the  spectral  estimate  than  Bartlett's  window. 

2.9.4  Parzen's  Lag  Window 

Definition  2.9.4  The  function  defined  by 


!-6[  jftt-  r * 6[  Ifti.  ]> . |u|  <| 


Vu)  *1  2d  - Y 3‘ 


< u < M 


, |u|  > M 


(2.9.8) 


is  called  Parzen's  lag  window.  The  Fourier  transform  of  equation  (2.9.8)  is 
called  Parzen's  spectral  window  and  is  given  by 


Wp(f)  « f M 


3 u rSin(fTTM/2)i» 


< f < 


(2.9.9) 
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The  above  window,  which  was  introduced  in  the  early  1960's,  gives  much  wider 
lobes  and  its  power  is  not  concentrated  primarily  at  lower  frequencies  as  in 
the  previous  window. 

2.10  Some  Remarks  Concerning  the  Lag  Windows 

Scientists  who  have  been  involved  in  choosing  the  shape  of  a lag  window, 
w(u),  have  taken  into  consideration  the  fact  that  the  spectral  window,  W(f) 
(that  is,  the  Fourier  transform  of  the  lag  window),  should  be  concentrated 
near  the  zero  frequency.  Blackman  and  Tukey,  [13],  looking  at  the  problem 
from  the  communications  engineering  point  of  view,  almost  identified  it  with 
that  of  choosing  the  intensity  distribution  along  an  antenna,  so  that  the 
variation  will  fall  in  a narrow  beam.  The  principal  maximum  and  the 
subsidiary  extrema  of  W(f)  are  called,  respectively,  main  and  side  lobes. 

A window  should  be  an  even  function  so  that  it  can  equally  treat  positive  and 
negative  values  of  the  spectral  density  function  on  both  sides  of  a given 
point  of  the  time  series.  It  should  integrate  to  unity,  that  is, 

0*(f)df  = 1 , (2.10.1) 

and  should  achieve  a maximum  value  at  the  frequency  f = 0.  That  is, 

|W(f) I < W(0),  for  all  f. 

It  should  be  concentrated  as  much  as  possible  about  f = 0 so  that  the 
behavior  of  the  spectral  density  function  is  concentrated  as  much  as  possible 
in  that. neighborhood. 

It  is  my  opinion  that  there  is  no  agreed  valid  criterion  for  comparing 
the  degree  of  concentration  of  any  window.  One  criteria  could  be  the  ratio 
of  the  size  of  the  second  largest  peak  to  the  size  of  the  largest  peak. 
However,  again  this  would  be  powerful  only  in  the  case  where  the  second 
largest  peak  would  occur  at  the  same  point.  This  fact  explains  why  one  has 
to  consider  all  the  different  windows,  in  addition  to  the  most  popular,  in 
one's  search  for  the  most  appropriate  case. 

For  the  main  lobe  of  W(f)  to  be  concentrated,  the  graph  of  w(u)  should 
be  flat  due  to  the  way  the  two  concepts  are  related.  Also,  for  the  side 
lobes  to  be  small,  w(u)  should  be  smooth  and  should  not  change  rapidly  as  in 
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the  case  of  the  rectangular  window.  Therefore,  one  should  compromise.  The 
above  authors'  investigations  have  been  done  along  the  lines  of  compromise, 
and  as  a result,  there  exist  numerous  windows  among  which  to  choose. 

Taking  Bartlett's  spectral  window,  Wg(f),  as  an  example,  we  find  that 

when  it  is  graphed  against  frequency,  it  is  found  to  be  symmetric  about  the 

1 2 3 

origin  and  has  zeros  at  f = ±rr,  ±g,  ±^,  ...  . 

The  distance  between  the  first  zeros  on  either  side  of  the  origin  is 
called  base  width.  The  base  width  for  Bartlett's  window  is  equal  to  ^ . It 
is  inversely  proportional  to  M and  also  to  the  variance.  On  the  other  hand,, 
by  increasing  the  base  width,  the  bias,  8(f),  increases  as  well.  Thus,  one 
is  forced  to  compromise  between  bias  and  variance  in  choosing  a particular 
window. 

The  rectangular  window  is  more  concentrated  about  the  center  frequency 
than  any  of  the  other  windows  under  consideration.  Nevertheless,  although 
it  has  the  smallest  bandwidth,  which  implies  small  bias,  it  also  has  the 
largest  side  lobes.  This  makes  it  very  impractical.  The  first  side  lobe  is 
about  1/5  of  the  height  of  the  main  lobe  which  shows  unrealistic  characteri- 
zation of  the  estimate  of  the  power  spectrum. 

Tukey  claims  that  the  window  he  proposed  with  Blackman,  the  use  of  which 
is  called  "Hanning"  after  the  Austrian  meteorologist  Ju.ius  Von  Han,  is 
simple  and  convenient.  Two  facts  about  it  are:  that  the  main  lobe  is  four 
times  as  wide  as  the  side  lobes;  and  that  the  side  lobes  are-  1%  or  2%  of  the 
height  of  the  main  lobe. 

The  above  remarks  will  be  extensively  examined  with  the  analysis  of  the 
ionospheric  information  in  Section  A. 
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3.  MODELING  AND  ANALYSIS  OF  IONOSPHERIC  INFORMATION 


3.1  INTRODUCTION 

HF  communications  are  provided  by  systems  not  specifically  limited  by 
1 ine-of-sight,  extended  distance,  or  intervening  terrain  obstacles.  However, 
ionospheric  disturbances,  both  natural  and  man-made,  complicate  the  use  of  the 
HF  media,  [2].  The  HF  communicator  had  as  his  only  propagation  aid,  the 
monthly  predictions  for  undisturbed  conditions  prepared  three  months  in  advance 
by  the  Department  of  Commerce,  and  distributed  by  the  U.  S.  Army  Strategic 
Conmuni cations  Command.  While  valid  for  long  range  planning,  they  do  not 
account  for  diurnal  variations  or  disturbed  ionospheric  conditions  that  may 
degrade  communications.  It  is,  therefore,  necessary  to  develop  a system  to 
provide  tactical  communicators  with  propagation  predictions  in  near  real-time, 
and  prepared  specifically  for  medium  and  long  range  distances. 

The  aim  of  this  section,  therefore,  is  two-fold: 

a.  to  introduce  a widely  accepted  statistical  concept  for  the  prediction 
of  oblique  incidence  soundings  including  application  for  the  prediction  of 
vertical  incidence  soundings,  and 

b.  to  develop  more  suitable  statistical  models  to  forecast  either  the 
oblique  or  vertical  incidence  soundings  over  specific  paths  or  at  specific 
terminals,  one,  two,  three,  ....  k time  slots  ahead,  beginning  with  a certain 
origin. 

It  is  shown  that  for  a 500  Km  path,  both  oblique  and  vertical  incidence 
recordings  are  non-stationary  stochastic  realizations.  That  is,  they  form  a 
discrete  time  series  that  is  not  in  statistical  equilibrium.  A procedure  is 
proposed  to  handle  this  type  of  information  and  to  investigate  the  possibility 
of  characterizing  the-  data  with  either  an  autoregressive  process,  a moving 
average  model,  or  a mixture  of  autoregressi ve-moving  average  processes. 

A systematic  presentation  of  recording  ionospheric  soundings  for  the 
purpose  of  forecasting  is  given  in  section  3.2.  An  autoregressive  model  has 
been  developed  for  the  discrete  realization  representing  the  13th  day  vertical 
incidence  (VI)  and  oblique  incidence  (01)  critical  frequencies  observed  for 
the  500  Km  path.  Fort  Monmouth,  N.  J.  - Fort  Drum,  N.  Y.,  in  section  3.3.  The 
complete  procedure  of  fitting  the  models  is  given,  along  with  the  associated 
confidence  Intervals.  Also  included  in  this  section  is  the  development  of 


47 


another  autoregressive  model  that  characterizes  the  behavior  of  the  18-day 
overall  VI  and  01  soundings  for  the  experiment.  A prediction  process  based 
on  the  widely  accepted  Ames-Egan  approach,  [1],  is  given  in  section  3.4.  In 
this  section,  the  limitations  of  the  process  and  suggested  improvements  are 
covered.  Finally,  in  section  3.5,  a summary  and  conclusion  are  presented. 

3.2  DESIGN  OF  THE  EXPERIMENT 

To  accomplish  the  task  of  developing  both  a functional  relationship 
between  01  and  VI  maximum  observed  frequencies  (MOF)  and  a forecasting  model, 
the  Communications/ADP  Laboratory,  of  the  U.  S.  Army  Electronics  Command,  had 
been  involved  in  an  extensive  collection  of  VI  and  short-path  01  ionospheric 
data  at  three  different  distances  with  Fort  Monmouth,  N.  J.,  as  the  base 
station.  Experimentation  was  performed  in  the  2-16  MHz  range,  using  two 
ionosondes,  one  as  a fixed  terminal  and  the  other  as  a mobile  terminal,  as 
shown  in  figure  3.1.  The  mobile  terminal  was  situated  at  Fort  Oix,  N.  J., 
establishing  a 50  i<m  path;  at  Aberdeen  Proving  Ground,  Md.,  to  establish  a 
200  Km  path;  and  at  Camp  Drum,  N.  Y.,  to  establish  a nominal  500  Km  path 
(figure  3.2).  This  section  will  address  only  the  500  Km  experiment. 

Each  terminal  made  scheduled  soundings  every  ten  minutes  for  an  18-day 
experiment.  While  the  fixed  terminal  was  transmitting  and  receiving  its  own 
signal,  the  mobile  terminal  would  simultaneously  receive  the  same  transmis- 
sions. The  same  procedure  was  followed  for  the  mobile  terminal  with  respect 
to  the  fixed  terminal  (figures  3.1  and  3.2).  Both  ionosondes  were  synchro- 
nized to  the  WWV  (HF)  time  standard  (National  Bureau  of  Standards)  so  that  the 
"remote"  sounder  scans  would  be  precise  with  the  Fort  Monmouth  terminal.  The 
number  of  days  each  experiment  was  performed  has  no  significance  with  respect 
to  the  results  obtained,  but  was  a matter  of  funding.  The  basic  instruments 
used  were  two  Granger  Associates  Model  3905-5  Ionospheric  Sounders,  matched 
with  wide  responde  delta  antennas. 

The  frequency  range  of  the  ionosondes  was  limited  from  2-16  MHz,  in  three 
octaves,  with  400  discrete  frequency  channels  per  octave.  Transmissions 
consisted  of  successively  "stepping"  through  the  channels  of  each  octave  with 
a pulse  width  of  100  micro-seconds  to  maximize  the  system  sensitivity.  The 
data  is  a recording  of  the  time  delay  from  ionosonde  to  ionospheric  reflecting 
layer  and  return.  Time  delay  is  a measure  of  the  virtual  height  of  reflection 
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from  the  layer.  The  trace  of  the  returned  pulse  on  a scale  of  frequency 
versus  time  delay  (virtual  height)  is  the  ionogram  record,  figure  3.3. 

Ionogram  records  of  the  data  were  taken  on  35im  film  at  Fort  Monmouth  and  on 
light  sensitive  oscillograph  paper  at  the  remote  terminals.  After  collection 
and  development,  the  ionograms  were  scaled  for  the  extraordinary  critical 
frequencies,  fxF ^ , as  shown  in  figure  3.4.  The  f Fg  data  was  then  compiled 
for  computer  analysis  and  for  comparison  between  the  observed  VI  and  observed 
01  critical  frequencies. 

The  experiment  results  were  dependent  upon  ionospheric  conditions  and 
man-made  noise.  Conditions  were  characterized  by  the  Space  Disturbance  Fore- 
cast Center,  Institute  of  Telecommunication  Sciences,  Boulder,  Colorado,  as 
generally  undisturbed,  but  some  interference  occurred.  Some  data  (ionograms) 
were  unreadable  due  to  man-made  noise,  and  solar  and  geomagnetic  activity. 

For  those  few  records  that  were  unreadable  (though  signal  was  detected), 
simulated  data  was  prepared.  The  occurrence  of  obscured  data  was  negligible 
over  the  experiment. 

3.3  FORECASTING  MODELS  FOR  IONOSPHERIC  SOUNDINGS 

In  this  section,  we  shall  propose  a step-by-step  procedure  in: 

a.  identifying  and  filtering  the  ionospheric  information, 

b.  fitting  the  most  appropriate  model  to  the  data, 

c.  applying  the  backwards  filter  technique  and  diagnostic  checking,  and 

d.  forecasting  and  updating  of  the  appropriate  model. 

Specific  reasons  and  mathematical  formulations  for  this  sequence  of  steps 
are  stated  in  Section  2.  Therefore,  the  basic  idea  and  philosophy  of  this 
section  is  the  implementation  of  the  procedural  approach  proposed  in  this 
study  for  analyzing  non-stationary  information.  The  approach  has  yielded 
better  results  for  fitting  and  forecasting  than  those  obtainable  with  other 
methods,  [1],  &4],  [15],  [16],  [17 3 » Cl8]»  that  exist  in  the  literature. 

In  the  previous  section,  the  manner  in  which  our  time  series  were 
generated  was  outlined  for  the  vertical  incidence  (VI)  and  oblique  incidence 
(01)  ionospheric  soundings.  Specifically,  we  have  available  18  days  of  data 
where  each  day  is  divided  into  144  time  slots.  The  data  was  taken  over 
consecutive  days  under  the  same  conditions;  however,  due  to  diurnal  ionospheri 
changes,  there  is  an  inherent  and  uncontrollable  diurnal  variation  in  the 
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soundings.  Thus,  the  observed  data  is  a stochastic  realization  that  consists 
of  144  data  points  with  the  corresponding  time  of  day.  In  view  of  the  diurnal 
differences  in  the  data,  it  was  .appropriate  to  choose  one  of  the  18  days  at 
random  for  each  of  the  observed  ionospheric  series  for  the  time  series 
analysis.  Furthermore,  we  felt  that  it  was  appropriate  to  also  consider  the 
average  of  the  18  days  for  the  time  series  analysis.  The  13th  day  (April  5, 
1971)  was  chosen  at  random  for  the  statistical  analysis.  This  amounts  to 
doing  a very  explicit  time  series  analysis  with  respect  to  identifying  the 
appropriate  stochastic  process  that  characterizes  these  four  phenomena.  We 
shall  begin  by  following  the  procedural  approach  detailed  in  Section  2,  which 
we  believe  is  a very  good  procedure  for  both  short-term  and  near-real -time 
forecasting  for  this  type  of  information. 

3.3.1  Model  Identification  and  Filtering 

The  initial  step  in  the  procedural  approach,  as  was  stated  previously, 
is  the  identification  of  the  stochastic  realization.  Primarily,  we  ask,  "is 
the  information  that  we  have  to  analyze  stationary  or  non-stationary?"  By 
stationary  we  mean  that  the  data  will  be  in  equilibrium  around  a constant 
mean  without  any  trends.  The  original  information  shown  by  figures  3.5,  3.6, 
3.7,  3.8,  which  represent,  respectively,  the  13th  day  observed  VI,  mean  VI, 
13^  day  observed  01,  and  mean  01  data,  visually  represent  non-stationary 
stochastic  realizations.  In  order  to  justify  this  fact,  the  sample  autocor- 
relation function  was  plotted  and  a statistical  test  was  performed  to  see  if 
there  were  any  non-linear  components  in  the  data.  Figures  3.9,  3.10,  3.11, 
and  3.12  show  that  in  all  cases  the  sample  autocorrelation  function  does  not 
dampen  out  very  rapidly.  This  is  certainly,  as  stated  in  section  2.4,  an 
indication  that  there  are  non-linear  components  within  the  observations. 
Secondly,  Kendall's  Tau  statistic,  [19]*  for  each  series  at  the  a = .05  level 
of  significance  (as  shown  in  table  3.1,  p.  63)  indicates  that  there  are  non- 
linear trends. 

* Chapter  5 
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3.5  OBSERVED  VERTICAL  INCIDENCE  CRITICAL  FREQUENCIES 
FOR  FORT  MONMOUTH,  N.  J.,  5 APRIL  1971 


FIGURE  3.7  OBSERVED  500  Km  OBLIQUE  INCIDENCE  CRITICAL  FREQUENCIES  BETWEEN 
FT.  MONMOUTH,  N.  J.,  AND  FT.  DRUM,  N.  Y.,  5 APRII  1971 


FIGURE 3. 9 Mean  500  Km  Oblique  Incidence  Critical  Frequencies  Between 

Ft,  Monmouth,  N.  J.,  and  Ft,  Drum,  N,  Y.,  22  March-10  April  1971 
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FIGURE  340  SAMPLE  AUTOCORRELATION  OF  THE  MEAN  VERTICAL  INCIDENCE  CRITICAL  FREQUENCIES 

FOR  FORT  MONMOUTH,  N.  J,,  22  MARCH  - 10  APRIL  1971 
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FIGURE  3.11  SAMPJ.E  AUTOCORRELATION  OF  THE  OBSERVED  500  Km  OBLIQUE  INCIDENCE  CRITICAL  FREQUENCIES 

BETWEEN  FT.  MONMOUTH,  N,  J.  AND  FT.  DRUM,  N.  Y.,  5 APRIL  1971 
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Table  3.1  Kendall's  Tau  Statistics  for  Trend 


Ionospheric  Series 

Calculated 

Statistic 

Z.  05 

Dec  i s i on 

13tl1  day  VI  observed 

5.192 

±1.645 

reject  HQ 

18-day  VI  averaged 

4.874 

±1.645 

reject  Hq 

13th  day  01  observed 

5.641 

±1.645 

reject  H 

0 

18-day  01  averaged 

4.964 

±1.645 

reject  Hq 

Thus,  we  are  certain  that  the  ionospheric  information  is  not  in  statis- 
tical equilibrium.  The  next  step  in  the  proposed  procedure  is  to  develop  a 
filter  that  will  eliminate  the  non-stationary  components  from  the  series.  The 
initial  step  towards  this  end  is  to  implement  a first  difference  filter, 
namely: 

yt  3 xt  - t 3 1 , ...,  144, 

to  each  of  the  four  realizations.  Upon  applying  this  filter,  the  sample  auto- 
correlation functions  were  computed  for  the  four  cases  and  are  shown  in 
figures  3.13,  3.14,  3.15,  and  3.16.  It  is  evident  by  inspecting  these  autocor- 
relation functions  that  the  13th  day  observed  01  and  VI  information  dampen  out 
fairly  rapidly  about  the  zero  axis  except  at  the  zero  point.  Kendall's  Tau 
test  was  conducted  on  the  filtered  data.  Results  indicate  that  in  both  of 
these  cases  the  non-stationary  components  had  been  removed  as  shown  in 
table  3.2. 


Table  3.2  Kendall's  Tau  Statistics  for  the 
First  Difference  Filtered  Data 


Ionospheric  Series 

Calculated 

Statistic 

Z.05 

Decision 

13^  day  VI  observed 

-0.379 

±1.645 

accept  Hq 

18-day  VI  averaged 

-4.577 

±1.645 

reject  H 

0 

13tfl  day  01  observed 

-0.574 

±1.645 

accept  H 

0 

18-day  01  averaged 

-4.616 

±1.645 

reject  Hq 
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Thus,  Kendall's  Tau  test  and  the  visual  interpretation  of  the  sample 

i,  u 

autocorrelation  function  indicate  that  we  have  reduced  the  two  13  day 
observed  stochastic  realizations  (01  and  VI)  into  stationary  series.  Thus, 
they  are  in  the  proper  form  to  proceed  in  identifying  the  appropriate  model 
that  characterizes  their  behavior.  However,  since  the  average  of  the  18  days 
of  01  and  VI  soundings  do  not  pass  the  criteria  for  Kendall's  Tau  test,  then 
one  concludes  that  trends  still  exist  in  these  series.  Furthermore,  the 
sample  autocorrelation  functions  in  these  cases  do  not  indicate  that  there 
is  rapid  enough  dampening  to  insure  stationarity  (see  figures  3.14  and  3.16). 
Therefore,  we  must  proceed  to  use  a second-difference  filter,  namely: 

wt  * xt  “ 2xt-1  + xt-2’  * = • 

As  shown  in  figures  3.17  and  3.18,  the  sample  autocorrelation  functions 
of  the  second-difference  averages,  the  dampening  behavior  seems  to  have  been 
improved.  Furthermore,  Kendall's  Tau  test  (as  shown  in  table  3.3  below)  seems 
to  have  reduced  the  original  non-stationary  realizations  into  proper  form; 
thus,  we  can  proceed  to  model  the  information. 


Table  3.3  Kendall's  Tau  Statistics  for  the 
Second  Difference  Filtered  Data 


Ionospheric  Series 

Calculated 

Statistic 

Z.  05 

Decision 

18-day  VI  averaged 

0.053 

±1.645 

accept  Hq 

18-day  01  averaged 

-0.093 

±1.645 

accept  H 

0 

It  is  appropriate  at  this  time  to  emphasize  that  the  13th  day  observed  01  and 
VI  data  needed  only  the  first-difference  filter.  The  non-stationarities  that 
occur  in  each  of  the  18  days  seemed  to  have  increased  the  non-stationary 
components  during  the  averaging  process,  requiring  a second-order  filter. 
However,  as  shown  in  figures  3.17  and  3.18,  the  improvement  in  the  sample 
autocorrelation  function  was  slight.  Since  we  feel  that  there  is  a stronger 
decision  with  the  second-difference  filter  for  the  averaged  series,  the 
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procedural  approach  will  be  continued  with  the  first-order  filter  in 
analyzing  the  13th  day  observed  data,  and  with  the  second-order  filter  in 
analyzing  the  averaged  information. 

3.3.2  Fitting  The  Model 

Now,  having  reduced  the  stochastic  realizations  into  stationary  form, 
the  next  step  in  the  procedural  approach  is  the  fitting  process.  Specifically, 
we  are  interested  in  identifying  the  appropriate  model,  and  secondly,  having 
identified  the  model  as  either  autoregressive  (AR),  moving-averages  (MA),  or 
a mixture  of  the  two  (ARMA),  the  order  of  the  particular  process  must  be 
determined.  In  deciding  on  the  order  of  a particular  process,  one  has  to  be 
concerned  with  estimating  the  parameters  that  are  associated  with  the  specific 
model.  Specifically,  with  regard  to  the  four  ionospheric  series,  we  have  a 
simultaneous  investigation  of  both  the  model  that  fits  that  data  and  the  esti- 
mation of  the  parameters  associated  with  the  model.  As  indicated  in  the 
recommended  procedure  for  identifying  the  model,  graphic  displays  have  been 
structured  as  shown  by  figures  3.19,  3.20,  3.21,  and  3.22,  for  each  of  the 
four  realizations  that  give  the  residual  variance  of  each  of  the  realizations 
as  a function  of  the  order  of  the  three  models  involved.  The  decision  as  to 
which  model  best  characterizes  these  series  resulted  from  the  criterion  of 
minimum  residual  variance  as  stated  in  section  2.4.2. 

As  shown  in  the  four  graphs,  the  order  (m,q)  refers  to  the  order  m of  the 
autoregressive  process  and  the  order  q of  the  moving  averages  process.  Thus, 
for  example,  (3,0)  is  a purely  AR  process  of  order  3;  (0,2)  is  a purely  MA 
process  of  order  2;  and  (3,2)  is  a mixture  of  a third  order  AR  and  a second 
order  MA  process.  As  depicted  in  the  figures,  with  respect  to  the  residual 
variance,  it  is  clear  that  the  "best"  models  that  should  be  selected  are: 

i)  (2,0)  process  for  the  13th  day  observed  VI  soundings, 

ii)  (3,0)  process  for  the  18-day  mean  VI  soundings, 

iii)  (3,0)  process  for  the  13tfl  day  observed  01  soundings, 

iv)  (3,0)  process  for  the  18-day  mean  01  soundings. 

It  will  be  shown  later,  that  the  difference  equations  that  have  been 
identified  above  are  the  most  appropriate  ones  to  characterize  our  data.  Of 
course,  in  reaching  the  decisions  with  respect  to  the  selection  of  these 
models,  the  principle  of  "parsimony"  was  also  taken  into  consideration. 
Specifically,  this  means  that  if  there  was  a lower  order  model  (refer  to 
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FIGURE  3.19  MODEL  ORDER  ys.  RESIDUAL  VARIANCE  FOR  THE  VERTICAL  INCIDENCE  IONOSPHERIC  DATA 

FT.  MONMOUTH,  N.  J.,  FOR  5 APRIL  1971 


Order  - m,q 

FIGURE  3.20  MODEL  ORDER  vs.  RESIDUAL  VARIANCE  FOP  THE  VERTICAL  INCIDENCE  IONOSPHERIC  DATA, 

FT.  MONMOUTH,  N.  J.  - MEAN  FOR  22  MARCH-10  APRIL  1971 
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! MODEL  ORDER  y§,  RESIDUAL  VARIANCE  FOR  THE  500  Km  OBLIQUE  INCIDENCE  IONOSPHERIC  DATA 
BETWEEN  FT.  MONMOUTH,  N.  J.,  AND  FT,  DRUM,  N.  Y.'f  - MEAN  FOR  22  MARCH-10  APRIL  1971 


figures  3.19  through  3.22)  that  had  an  equal  or  slightly  higher  residual  vari- 
ance than  a higher  order  model,  the  lowest  order  model  would  be  selected 

because  the  number  of  parameters  and  convenience  that  such  a model  offers 

computationally  is  desirable,  and  will  usually  have  little  effect  on  the  near- 
real  -time  forecasts  of  the  data. 

The  parameter  estimates  of  the  specified  model  are  simultaneously 
obtained  with  the  determination  of  the  order.  To  calculate  the  least  squares 
estimates  of  the  parameters,  the  procedure  discussed  in  section  2.4  was  fol- 
lowed. The  estimates  of  the  AR  parameters,  a , and  those  of  the 

MA  process,  8-j , 82*  . ..,  8 > must  satisfy  the  stationarity  and  invertibi  1 ity 

properties,  respecti vely.  Briefly  stated,  for  each  model  (m,q)  considered 
(see  figures  3.19  through  3.22),  one  finds  the  stationary  region  for  the  m 
autoregressive  parameters  and  the  invertibi! ity  region  for  the  q moving- 
average  parameters,  and  forms  the  joint  region  for  both.  This  region  is  then 
"gridded"  over  all  parameters  and  the  residual  variance  is  computed  for  each 
point  (a.j , a^,  ...»  a^,  8^,  82,  ...,  6q).  One  then  finds  the  point  of  minimum 
residua!  variance.  Thus,  we  choose  the  model  (m,q)  and  the  parameters  associ- 
ated with  it,  which  result  in  minimum  residual  variance.  To  do  this  we  used  a 
grid  program  which  computed  the  residual  variance  for  all  possible  combina- 
tions of  parameters  for  each  of  the  three  processes.  Therefore,  for  the 
appropriately  filtered  series,  the  estimates  of  the  true  states  of  nature  are 
shown  in  table  3.4. 


Table  3.4  Approximate  Least  Squares  Estimates 
of  the  Best  Model  Parameters 


Model 

Series 

Order  (m,q) 

" “a- 

u 

al 

* 

a2 

a3 

i . 

13tl1  day  VI  observed 

(2,0) 

0.0020 

0.289 

0.395 

— 

ii . 

18-day  VI  averaged 

(3,0) 

0.0004 

-0.676 

-0.559 

-0.258 

iii . 

13tfl  day  01  observed 

(3,0) 

0.0000 

0.484 

-0.102 

0.266 

iv. 

18-day  01  averaged 

(3,0) 

0.0008 

-0.574 

-0.352 

-0.070 

These  results  yielded,  respectively,  the  following  difference  equations  in 
terms  of  the  first  difference  filtered  series,  yt>  and  the  second  difference 
filtered  series  w • 
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(3.3.1) 


i.  (yt  - .002)  = .289  (yt_]  - .002)  + .395  (yt-2  - .002)  + Zt 

ii.  (wt  - .0004)  = -.676.  (wt_1  - .0004)  - .559  (wt_2  - .0004) 

- .258  (wt_3  - .0004)  + Zt  (3.3.2) 

iii .  (yt  - 0.000)  = .484  (yt-1  - 0.000)  - .102  (yt_2  - 0.0000) 

-t-  .266  (yt_3  - 0.000)  + Zt  (3.3.3) 

iv.  (wt  - .0008)  = -.574  (wt_,  - .0008)  - .352  (wt_2  - .0008) 

- .070  (wt_3  - .0008)  + Zt  (3.3.4) 

3.3.3  Inserting  the  Backwards  Filter  and  Diagnostic  Check  of  the  Models 

The  next  step  in  the  procedural  approach  to  forecasting  is  the  implemen- 
tation of  the  backwards  filter  and  diagnostic  check  of  the  models.  Having 
selected  the  appropriate  stationary  stochastic  model  and  its  order,  a diag- 
nostic check  must  be  performed  to  determine  the  adequacy  of  the  models.  As 
indicated  in  Section  2,  if  the  original  information  was  filtered  to  put  it 
into  the  proper  form  to  perform  time-series  analysis,  we  must,  at  this  point 
incorporate  back  into  the  model  the  non-stationarities  that  the  filter  has 
eliminated.  That  is*  we  must  introduce  the  concept  of  the  backward  filter 
into  our  model.  This  backward  filtering  concept  is  very  important  because 
it  puts  back  into  the  model  some  of  the  basic  characteristics  that  the 
initial  data  contained  so  that  the  final  interpretation  of  the  observed 
realizations  would  be  more  meaningful.  Thus,  we  insert  the  appropriate 
filter,  either: 

*t  = xt  * xt-l  ’ 
or 

wt  * xt  - 2xt-l  + xt-2 

into  equations  (3.3.1),  (3.3.2),  (3.3.3),  and  (3.3.4). 

Given  below  are  the  forecasting  models,  having  been  modified  to  include 
the  filtering  concept.  That  is, 
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A 

i.  x^  = 1.289x^._i  + 0.106x.£_2  “ 0.395Xj._2  + 0.007  + (3.3.5) 

A 

ii.  xt  = 1.324xt_-j  - 0.207xt_2  + 0.184xt_3  - Q.043xt_4  - 0.258xt_5 

+ .001  + Zt  (3.3.6) 

/v 

iii.  x^.  = 1 .484x^._i  - 0.586x^  ^ + 0.368x^._2  - 0.266x^._4  + Z^  (3.3.7) 

iv.  xt  = 1 . 426xtl  + 0.204xt_2  - 0.060xt_3  - 0.212xt_4  - 0.070xt_5 

+ .0016  + Zt  (3.3.8) 


Having  formulated  the  above  difference  equations  for  the  stochastic  realiza- 
tions, the  next  step  is  to  investigate  the  goodness-of-fit  of  the  structured 
models.  This  is  done  as  outlined  in  section  2.4  by  calculating  the  residuals 
incurred,  that  is,  by  subtracting  the  modeled  series  from  the  observed  series. 

A 

In  other  words,  if  x^.  is  the  observed  series  and  x^  is  the  modeled  series, 
then  their  difference  r.,  is  the  residual  (equation  2.4.16);  that  is: 

/v  U 

rt  = xt  - xt*  The  residuals  should  behave  as  a purely  random  process,  with 
a zero  mean  and  a variance  in  the  order  of  1/n.  As  a first  step  in  simulat- 
ing the  x^.,  the  unknown  value  of  Zt  is  set  to  its  unconditional  expectation 

of  zero  in  equations  (3.3.5)  through  (3.3.8),  and  it  is  assumed  that  the 
values  of  xt_-j , xt_3>  • ••»  xt-m’  are  known.  As  shown  by  figures  3.23? 

3.24,  3.25,  and  3.26,  we  have  an  excellent  fit  of  the  estimated  models  with 

respect  to  the  observed  series.  One,  of  course,  can  perform  a statistical 
test  to  justify  the  fit.  That  is,  the  resulting  residuals  should  behave 
approximately  like  a purely  random  process;  in  other  words,  they  should  be 
normally  distributed  with  a mean  of  zero  and  a variance  of  1/n.  Clearly,  for 
this  sample  size,  the  variance  becomes  very  small.  Thus,  the  standard  devia- 
tion of  the  sample  autocorrelation  is  l//rT  = 1//144  = .0833.  The  95% 
confidence  intervals  of  the  sample  autocorrelation  r22(k)  are  ±1.96  (.0833) 

= ±0.163.  At  the  5%  level  of  significance,  one  could  expect  (.05)144  or  8 
out  of  the  sample  autocorrelations  to  lie  outside  the  confidence  interval. 

Only  two  autocorrelations  for  the  13th  day  observed  VI  and  01,  and 
two  autocorrelations  for  the  18-day  averaged  information  lie  outside 
of  the  confidence  interval  (see  tables  3.5,  3.6,  3.7  and  3.8).  Hence,  one 

*Note:  In  the  simulation  graphs,  the  lines  connecting  the  points  are  aids 
to  see  the  relationship  of  the  points. 
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FIGURE  3.23  SIMULATED  VI  SERIES  USING  THE  AUTOREGRESSIVE  MODEL  vs.  THE  OBSERVED  VI 
CRITICAL  FREQUENCIES  FOR  FT.  MONMOUTH,  N.J.,  5 APRIL  1971 


SIMULATED  500Km  01  SERIES  USING  THE  AUTOREGRESSIVE  MODEL  vs.  THE  OBSERVED 
500 Km  01  CRITICAL  FREQUENCIES  BETWEEN  FT.  MONMOUTH,  N.J.,  AND 
FT.  DRUM,  N.Y . , 5 APRIL  1971 
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Table  3.5  Sample  Autocorrelations  r (k)  for  the  Simulated  13th  Day 
Observed  VI  Data  With  a 95%  Confidence  Interval  of  ±0.163 


Lag  k 
(Time 
Slots) 

Sample 

Autocorrelation  rzz 

(k)  xl(f 

2 

1-10 

-14.1 

-31 .8 

-.003  - 

.001 

-.051 

-.040 

-.023  - 

.023 

-.037 

-.039 

11-20 

-.005 

-.033 

-.021 

.008 

.011 

-.036 

-.010 

.010 

-.038 

-.019 

21-30 

-.008 

-.026 

.007  - 

.008 

-.002 

-.008 

-.093  - 

.069 

.157 

-.096 

31-40 

-.097 

-.031 

.042  - 

.002 

.056 

.103 

.166 

.094 

.069 

.103 

41-50 

.072 

.045 

-.051 

.003 

.058 

-.033 

-.055  - 

.077 

-.089 

.155 

51-60 

.157 

.025 

-.010 

.017 

.039 

-.037 

-.027  - 

.020 

-.084 

-.111 

61-70 

-.120 

.095 

.122 

.029 

.003 

.072 

-.021 

.036 

-.004 

-.024 

71-80 

-.111 

-.153 

.002 

.042 

.039 

.055 

.026 

.021 

.025 

-.007 

81-90 

.009 

.005 

.038  - 

.023 

-.013 

.021 

.019 

.009 

-.099 

-.119 

91-100 

-.059 

.025 

-.018 

.035 

.057 

.023 

.009  - 

.009 

-.019 

-.074 

101-110 

.066 

.000 

-.039  - 

.035 

-.079 

-.067 

.041 

.077 

.031 

-.032 

m-120 

.019 

.073 

-.121  - 

.037 

-.056 

-.040 

-.097  - 

.119 

-.014 

.056 

121-130 

-.032 

-.087 

-.011  - 

.006 

-.023 

.018 

.041  - 

.096 

-.106 

-.041 

131-140 

-.038 

-.027 

.009  - 

.040 

.005 

-.051 

-.087  - 

.029 

.025 

-.036 

141-143 

-.455 

-.342 

1.147 
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Table  3.6  Sample  Autocorrelations,  r2Z ( k ) f°r  the  Simulated  18-0ay 
Averaged  VI  Data  With  a 95%  Confidence  Interval  of  ±0.163 


Lag  k 
(Time 
Slots) 

Sample 

Autocorrelation  r zz 

(k)  xlO* 

2 

1-10 

-13.5 

.84 

-17.1  - 

20.2 

-.050 

.032 

.067  - 

.005 

.035 

-.064 

11-20 

-.014 

-.011 

-.003 

.017 

.032 

-0.21 

.017  - 

.044 

-.002 

.048 

21-30 

-.037 

.028 

-.097 

.072 

-.068 

.180 

-.095  - 

.031 

-.078 

-.023 

31-40 

.060 

.046 

-.106 

.067 

.014 

.099 

.008  - 

.028 

.018 

-.059 

41-50 

.033 

-.006 

.008 

.005 

-.033 

.000 

-.019 

.031 

-.011 

.013 

51-60 

.000 

-.015 

.035  - 

.032 

.027 

-.043 

.004  - 

.013 

.013 

.012 

61-70 

-.016 

.001 

-.003 

.002 

-.009 

.033 

.015 

.016 

-.023 

-.029 

71-80 

-.059 

.058 

-.031 

.075 

-.003 

-.018 

.017  - 

.029 

-.011 

.019 

81-90 

1 -.011 

1 

.007 

-.011 

.012 

.020 

.005 

-.001  - 

.045 

.001 

-.013 

91-100 

.021 

-.004 

-.014  - 

.006 

.008 

.018 

.012  - 

.008 

-.009 

.011 

101-110 

! -.019 

.025 

-.037 

.002 

.018 

.013 

.003 

.013 

-.017 

.029 

111-120 

-.001 

-.004 

.000  - 

.018 

-.003 

-.016 

-.013 

.000 

-.005 

.015 

121-130  j 

i -.031 

.035 

-.022 

.031 

.014 

-.016 

-.009  - 

.021 

-.027 

.006 

131-140 

.006 

-.006 

.034  - 

.010 

.018 

.000 

-.014 

.098 

.028 

.028 

141-143 

.086 

-.417 

.206 

84 


Table  3.8  Sample  Autocorrelations,  rzz(k)  for  the  Simulated  18-Day 
Averaged  01  Data  With  a 95%  Confidence  Interval  of  ±0.163 


Lag  k 
(Time 
Slots) 


Sample  Autocorrelation  rzz(k)  xlO*-2 


1-10 

-18.9 

-6.59 

-19.1 

-5.33 

-.012 

.018 

.006 

.009 

-.013 

-.010 

11-20 

-.007 

.000 

.029 

-.016 

-.009 

-.005 

-.004 

.008 

-.010 

.019 

21-30 

-.005 

.006 

.000 

-.028 

.041 

-.017 

.011 

-.063 

-.001 

.004 

31-40 

.011 

.001 

-.025 

.024 

.015 

.036 

.009 

.000 

-.021 

-.022 

41-50 

.023 

.008 

.003 

-.019 

-.001 

-.015 

.000 

.000 

.012 

.008 

51-60 

.003 

-.014 

.000 

.011 

.000 

-.012 

-.018 

.014 

.000 

.007 

61-70 

-.019 

.004 

-.008 

.013 

.005 

-.004 

.008 

-.020 

.016 

-.018 

71-80 

-.008 

.013 

.019 

-.003 

.000 

-.012 

.007 

-.006 

.003 

-.006 

81-90 

.004 

-.001 

.005 

.002 

.006 

-.009 

-.002 

-.011 

.004 

.000 

91-100 

-.002 

-.007 

-.007 

.003 

.012 

.000 

.013 

-.012 

-.001 

-.002 

101-110 

-.004 

.001 

-.009 

-.018 

.022 

.004 

.015 

.018 

-.014 

.001 

111-120 

-.010 

.007 

-.004 

-.009 

-.005 

.003 

-.013 

.013 

-.009 

.011 

121-130 

-.015 

.008 

-.004 

.020 

-.005 

-.006 

-.009 

-.011 

-.009 

.011 

131-140 

-.004 

.017 

-.002 

-.007 

-.005 

-.006 

.012 

-.073 

-.272 

-.135 

141-143 

-.385 

1.17 

-.319 
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can  conclude  that  the  residuals  constitute  a purely  random  process  and  that 
the  fitted  models  give  satisfactory  representation  of  the  observed  series. 

There  is  a further  statistical  test  described  in  section  2.4.3  that 
reflects  on  the  accuracy  of  our  models.  From  equation  (2.4.17),  the  calcu- 
lated value  of  Q 'v  x£_m_q  for  the  first  K autocorrelations,  where  K = n/10, 
m = the  order  of  the  AR  process,  q - the  order  of  the  MA  process.  When  this 
value  is  tested  against  the  appropriate  x2  value,  and  if 

« < Xk-m-q  <>  - ’ 

then  the  fitted  model  is  adequate,  [12].  For  the  stationary  stochastic 
ionospheric  realizations, 

i.  for  the  13th  day  observed  VI  data,  Q = 1 6 . 67 < x^  (.975)  = 24.7 

ii.  for  the  18-day  averaged  VI  data,  Q = 1 2 . 72 < (-975)  = 24.7 

iii.  for  the  13th  day  observed  01  data,  Q = 21.29 < x|3  (.975)  =23.3 

iv.  for  the  18-day  averaged  01  data,  Q = 10.41  < x^  (-975)  = 24.7. 

Clearly,  then,  the  four  estimated  models  are  adequate. 

By  the  mean  approaching  zero,  it  is  meant  that  in  using  these  types  of 
models  there  is  a tendency  to  either  under-  or  over-forecast  a particular 
estimate  at  a given  time  slot.  However,  in  the  long  run,  if  the  over-  and 
under-estimates  are  averaged,  this  average  would  be  zero.  The  basic  idea  of 
the  "purely  random  process"  means  that  as  n + «,  the  variance  approaches 
zero.  This  is  not  the  case,  however,  because  if  the  variance  is  zero,  we 
would  have  a degenerate  phenomenon  where  the  mass  would  be  concentrated  at 
a point,  and  interpretation  would  be  impossible. 

3.3.4  Forecasting  and  Updating 

The  final  step  in  the  procedural  approach  to  time  series  modeling  is  to 
state  the  model  in  such  a form  that  forecasts  l steps  ahead  are  possible, 
where  £=1,2,  ....  n.  Given  below  are  the  four  estimated  ionospheric 
models  in  the  appropriate  form  for  i step-ahead  forecasts: 

i.  for  the  13th  day  observed  VI  Information,  xtU)  = 1.289xt  .j 

+ 0.106  tn_2  - 0.395  tn_3  ♦ .007  ♦ ZtH  (3.3.9) 


4 


A 

ii.  for  the  18  - day  averaged  VI  information,  xt^  3 1.324xt+A_j 
- 0.207xt+J_2  + 0.184Xt+g_3  - 0..043xt+)l_4  - 0.258xt+Jl_5 
+ .001  + Ztn  (3.3.10) 

1*h  A 

iii.  for  the  13  day  observed  01  information,  xtuj  * 

0.586xt+ft_^  + 0.368xt+)l_3  - 0.266xt+^4  + Zt+|l  (3.3.11) 

iv.  for  theT8  - day  averaged  01  information,  x^4j  3 1.426xt+jl_( 
0.204xt+jl_2  + 0.060xt+jl_3  - 0L212xt+jl_4  - 0.070xt+s_,j 
+ Ztn  + .0016  . (3.3.12) 


Of  course,  the  accuracy  of  these  models  will  be  much  better  for  small  £.  As 
l becomes  large,  i.e.,  Z » m + d + q,  where  d is  the  order  of  the  filter  and 
m and  q are  as  previously  defined,  the  accuracy  decreases  substantially. 
However,  one  can  forecast  any  number  of  steps  ahead  and  update  the  forecasts 
as  additional  information  becomes  available. 

To  illustrate  how  one  may  update  the  forecasts  for  a time  t (origin), 
suppose  that  a new  piece  of  data,  xt+1 , becomes  available.  With  the  new 
origin  at  time  t+1 , we  update  the  time  t forecasts  by  equation  (2.4.42)  as: 

A A, 

*t+l  = x^(4+l ) +■  , £=1 , 2,  ...,  11, 

A 

where  Zt+1  = xt+1  - x^(l)  and  9^  is  as  described  in  section  2.4. 

Given  in  tables  3.9,  3.10,  3.11,  and  3.12  below  are  Z steps  ahead  fore- 


casts (up  to  £=11)  at  an  arbitrary  t = 72  origin,  with  updating,  along  with 


Table  3.9  Forecasted  Values  of  the  13th  Day  Observed  VI  Series  at 
Origin  t = 72  and  Updating  Under  the  Assumption  x73  Becomes  Available 


TIME 

ACTUAL 

VALUE 

mmgmm 

be  a 

FORECAST 

95%  PROBABILITY 
LIMITS 

UPDATED 

FORECAST 

1150 

9.00 

— 

.... 

8.80 

i 

8.875 

± .404 

1210 

8.70 

2 

8.885 

± .596 

8.788 

1220 

9.00 

3 

9.022 

t .778 

8.889 

1230 

9.20 

4 

9.047 

± .962 

9.198 

1240 

9.40 

5 

9.064 

±1.139 

9.235 

1250 

9.60 

6 

9.086 

±1.310 

9.270 

1300 

9.70 

7 

9.106 

±1 .474 

9.302 

1310 

9.80 

8 

9.127 

±1.631 

9.331 

1320 

9.90 

9 

9.148 

±1.782 

9.359 

1330 

9.90 

10 

9.169 

±1.926 

9.386 

1340 

10.00 

11 

9.190 

±1.926 

9.411 

Table  3.10 

Forecasted  Values  of  the 

18-Day  Averaged  VI  Series  at 

Origin 

t * 72  and  Updating 

Under  the  Assumption  x73  Becomes 

Available 

ACTUAL 

LEAD 

95%  PROBABILITY 

UPDATED 

TIME 

VALUE 

TIME 

FORECAST 

LIMITS 

FORECAST 

1150 

8.79 

-- 

8.71 

1 

8.721 

± .214 

— 

1210 

8.90 

2 

8.729 

± .292 

8.879 

1220 

8.96 

3 

8.752 

± .371 

8.969 

1230 

9.00 

4 

8.749 

± .473 

8.949 

1240 

9.13 

5 

8.737 

± .562 

8.987 

1250 

9.16 

6 

8.745 

± .616 

9.005 

1300 

9.17 

7 

8.752 

± .641 

8.972 

1310 

9.19 

8 

8.752 

± .644 

8.902 

1320 

9.19 

9 

8.753 

± .656 

8.803 

1330 

9.24 

10 

8.758 

± .735 

8.868 

1340 

9.25 

11 

8.763 

± .735 

9.043 
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Table  3.11  Forecasted  Values  of  the  13th  Day  Observed  01  Series  at 
Origin  t = 72  and  Updating  Under  the  Assumption  Becomes  Available 


TIME 

ACTUAL 

VALUE 

mm 

FORECAST 

95%  PROBABILITY 
LIMITS 

UPDATED 

FORECAST 

1150 

9.20 

-- 

1200 

9.30 

i 

9.305 

± .586 

1210 

9.60 

2 

9.439 

± .791 

9.509 

1220 

9.60 

3 

9.467 

± 1.006 

9.547 

1230 

9.90 

4 

9.495 

± 1.229 

9.585 

1240 

10.10 

5 

9.541 

± 1.438 

9.651 

1250 

10.20 

6 

9.568 

£ 1.638 

9.678 

1300 

10.20 

7 

9.584 

± 1.831 

9.704 

1310 

10.20 

8 

9.601 

± 2.016 

9.721 

1320 

10.60 

9 

9.615 

± 2.193 

9.745 

1330 

10.70 

10 

9.624 

± 2.362 

9.754 

1340 

10.60 

11 

9.631 

± 2.362 

9.761 

Table  3.12 

Forecasted  Values  of  the  18-Day  Averaged  01  Series  at 

Origin 

t = 72  and  Updating  Under  the  Assumption  x^3  Focomes 

Available 

ACTUAL 

95%  PROBABILITY 

UPDATED 

TIME 

VALUE 

mp^i 

FORECAST 

LIMITS 

FORECAST 

1150 

9.26 

.... 

1200 

9.25 

i 

9.249 

£ .198 

1210 

9.42 

2 

9.260 

£ .287 

9.360 

1220 

9.56 

3 

9.274 

£ .395 

9.474 

1230 

9.55 

4 

9.283 

± .514 

9.483 

1240 

9.62 

5 

9.294 

£ .652 

9.594 

1250 

9.66 

6 

9.307 

£ .801 

9.347 

1300 

9.69 

7 

9.320 

£ .959 

9.360 

1310 

9.71 

8 

9.334 

£1.129 

9.384 

1320 

9.75 

9 

9.349 

£1.649 

9.399 

1330 

9.75 

10 

9.365 

£1.495 

9.425 

1340 

9.78 

11 

9.382 

£1.495 

9.442 
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It  is  clear  from  tables  3.9  and  3.11  that,  for  a small  number  of  steps,  i.e., 
I » 1,  2,  3,  our  forecasts  of  the  observed  series  are  quite  good.  For 
larger  L,  the  forecasted  values  generally  tend  to  underestimate  the  future 
values.  Updating  the  forecasts  in  these  cases  improves  the  estimates.  The 
l step  ahead  forecasts  for  the  18-day  averaged  series  are  also  very  good  for 
l 3 1,  2,  3,  steps  ahead.  The  necersity  for  updating  in  these  cases  is 
obvious  since,  here  again,  updating  yielded  an  excellent  gain  of  accuracy. 

3.4-  A Prediction  Process  Based  on  The  Sample  Autocorrelation  Function 

Ames  and  Egan,  [1],  have  developed  a prediction  process  based  upon  the 
sample  autocorrelation  function  to  make  short-term  predictions  for  the  maxi- 
mum observed  frequency  (MOF),  the  lowest  observed  frequency  (LOF),  and  other 
high  frequency  (HF)  propagation  parameters  dealing  with  the  ionosphere.  The 
acquisition  of  the  data  for  their  model  was  similar  to  the  500  Km  experiment 
described  in  section  3.2. 

Their  data  was  collected  from  oblique  incidence  soundings  which  were 
transmitted  every  10  minutes  from  the  U.  S.  Naval  Station  at  Lualualei, 
Hawaii,  and  received  at  Palo  Alto,  California,  a distance  of  about  4400  Km. 
The  analysis  was  based  on  ionograms  recorded  from  January  7 through  March  12, 
1965.  For  each  24-hour  period,  there  were  144  increments  and,  in  duration, 
there  were  45  days  of  data  available. 

For  each  10-minute  increment  of  the  24-hour  period,  the  mean  and  the 
standard  deviation  were  calculated  for  the  MOF  and  LOF.  It  was  pointed  out 
that  the  standard  deviation  represents  the  variation  of  the  MOF  and  LOF  from 
day  to  day  within  each  particular  time  slot  and  not  fluctuations  with  time 
on  a single  day.  It  was  discovered  that  the  standard  deviation  of  the  MOF 
was  higher  during  the  day  than  at  night  by  a factor  of  2.17,  while  the  MOF 
ratio  between  these  periods  was  2.4  to  1.  The  plot  of  the  MOF  vs.  standard 
deviation  further  demonstrated  the  dependence  of  the  standard  deviation  upon 
the  MOF. 

Next,  the  sample  autocorrelation  function  of  the  MOF  was  calculated  with 
values  from  five  previous  10-minute  periods  and  then  plotted  as  hourly 
averages.  It  was  found  that  the  autocorrelation  of  the  MOF  was  high  during 
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the  morning  and  evening  ionospheric  transition  periods.  The  autocorrelation 
is  a measure  of  the  current  stability  relative  to  the  mean-value  function, 
whereas  the  standard  deviation  is  a measure  of  the  day-to-day  variance. 

3.4.1  Short-Term  Pr  jiction 

Ames  and  Egan  fitted  forecasting  models  to  the  maximum  observed 
frequency  series  and  to  the  lowest  observed  frequency  series  to  be  able  to 
predict  eight  successive  10-minute  periods  following  each  observation.  Their 
prediction  (or  forecasting)  model  is,  [1], 

x(t+T)  = x(t  +■  T)  + [x(t)  - x(t)]  p(t,T)  , (3.4.1) 

where  t is  the  present  time;  T is  the  lead  time;  x(t)  is  the  observed  value 
for  time  t;  x(t)  is  the  mean  value  for  the  time  slot  t;  x(t  + T)  is  the  long- 
term mean  value  for  a prediction  T minutes  in  the  future;  p(t,T)  is  the 

/V 

sample  autocorrelation  between  values  at  t and  t + T;  and  x(t+T)  is  the 
future  value  of  the  MOF  or  the  LOF  predicted  for  any  of  eight  10-minute 
periods  following  each  observation.  That  is,  equation  (3.4.1)  states  that 
a prediction  for  T minutes  in  the  future  consists  of  the  long-term  mean  value 
for  that  time  plus  a weighted  term  correcting  it  for  the  present  difference 
between  observed  and  average  values. 

The  input  data  Ames  and  Egan  used  in  equation  (3.4.1)  consisted  of 
current  MOF  and  LOF  running  averages  and  standard  deviations  determined  only 
from  previous  data  (except  for  the  first  few  days  of  "start  up"),  and  values 
of  autocorrelation  derived  from  all  the  data.  The  MOF  and  LOF  running 
averages  were  computed  and  stored  separately  for  each  of  the  10-minute  incre- 
ments of  the  24-hour  day;  the  standard  deviations  were  computed  for  each 
10-minute  increment  and  then  combined  into  hourly  averages.  These  values 
were  then  updated  when  a new  measurement  was  realized  by  the  running  average 
function: 

xn  = ^ ^ Vl  + T = ^ ^ X1  + ^ 1 Xi  * 

(3.4.2) 

It  was  stated  that  for  the  means,  xn  represented  an  observed  value  and  for 
the  standard  deviation,  it  represented  the  absolute  value  of  the  difference 
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between  an  observed  value  and  the  corresponding  mean.  The  value  of  K was 
set  equal  to  10  for  the  means,  and  100  for  the  standard  deviations. 

Data  was  used  from  the  full  experimental  period  in  calculating  the 
sample  autocorrelation  function  in  order  to  obtain  relatively  smooth  esti- 
mates of  the  autocorrelation.  The  following  equation  was  used  to  estimate 
the  autocorrelation: 


p(t,T) 


2N  [x. (t)  - x(t)][x.(t  - T)  - x{t  - T)] 

, Izl 1 I 

AN  [x.(t)  - x(t)]2  IN  [ x . ( t - T)  - x(t  - T) ]: 
/ i»l  i i=l  1 


(3.4.3) 

where  N is  the  number  of  days  measurements  were  available.  In  real-time, 
values  of  the  autocorrelation  could  be  computed  once  each  month  and  then  when 
new  values  are  realized,  they  could  be  updated  using  equation  (3.4.2). 

The  24  hourly  values  of  the  autocorrelation,  p,  for  a delay  of  10 
minutes  plus  a set  of  corresponding  time  constants,  t,  derived  from  the  least 
squares  fit  to  the  p data  between  delays  of  10  to  60  minutes,  were  used  as 
input  to  the  prediction  process.  For  delays  greater  than  10  minutes,  values 
of  p were  calculated  from: 


p = p1q  exp-(T  - 1Q)/t 


(3.4.4) 


The  values  of  p calculated  from  equation  (3.4.4)  were  found  to  be  a few 
percent  larger  than  the  observed  values  of  autocorrelation.  Ames  and  Egan 
concluded  that  the  smoothing  of  the  measured  p values  over  1-hour  periods 
plus  the  abstracting  of  all  delays  greater  than  10  minutes  into  corresponding 
time  constants  appears  to  have  reduced  the  unrealistic  benefit  from  this 
partial  view  of  the  future  to  a negligible  amount. 

Ames  and  Egan  express  the  expected  error  in  predicting  with  the  process 
(3.4.1)  in  terms  of  the  autocorrelation.  The  expected  error  is  given  by: 


ap(t  + T) 

of t + T) 


(t,  T) 


(3.4.5) 


The  above  equation  predicts  the  degree  to  which  the  variance  of  the  observed 
values  about  the  corresponding  predicted  values  will  be  less  than  that  of  the 
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observed  values  about  the  mean.  It  was  concluded  that  to  make  a substantial 
reduction  in  ap  relative  to  a,  say  one-half,  p must  be  greater  than  0.86. 
However,  it  was  stated  that  p. edictions  should  be  considered  "useful"  for 
values  of  p down  to  about  0.5  since  even  though  a is  not  much  reduced,  the 
predicted  value  itself  is  adjusted  from  the  long-term  mean  by  one-half  the 
presently  observed  difference. 

3.4.2  Updating  with  Fourier  Coefficients 

In  a later  article,  Ames,  Egan,  and  MacGinitie,  [20]>  use  Fourier 
coefficients  to  update  the  data  to  eliminate  the  random  variability  in  the 
input  data.  The  random  variability  caused  unavoidable  growth  of  irregular- 
ities in  the  diurnal  (daily)  curves  of  the  long-term  function.  Thus,  instead 
of  using  the  technique  of  running  averages  (equation  3.4.2)  to  update  the 
hourly  averages  when  a new  measurement  is  realized,  the  long-term  data  is 
converted  to  a limited  number  of  Fourier  coefficients  from  which  the  desired 
values  are  found  as  needed.  For  the  mean  values,  eight  harmonics  were  used; 
for  the  standard  deviations,  four  harmonics  were  used. 

The  following  equation  was  used  to  derive  the  Fourier  coefficients  from 
the  set  of  data  with  a 24-hour  period  of  144  10-minute  ir> tervals: 

an(NEW)  = an(0LD)  + e f [x(t)  - x(t)]  cos  (^)  , (3.4.6) 

where  e = 0.0944  and  T = 144.  The  value  of  e was  chosen  so  as  to  allow 
the  Fourier  coefficient  to  follow  long-term  ionospheric  changes  with  a time 
constant  of  approximately  11  days. 

Ames,  Egan,  and  MacGinitie  comment  that  this  conversion  to  Fourier 
coefficients  not  only  improves  the  prediction  quality,  but  also  substantially 
reduces  the  required  amount  of  computer  memory  capacity.  Also,  the  data 
storage  requirement  is  reduced  by  the  approximation  of  the  autocorrelation 
function  by  a decaying  exponential. 

3.4.3  The  Autocorrelation  Function 

Ames  and  Egan  throughout  their  paper  take  the  position  that  the  auto- 
correlation that  they  have  calculated  is  the  true  state  of  nature,  that  is, 
that  they  have  the  true  parameter  value,  not  just  an  estimate  of  it.  This 


is  most  clearly  shown  in  their  constant  reference  to  the  calculated 
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autocorrelation  as  p,  not  p,  or  r,  which  has  become  the  widely-accepted 
symbol  for  the  sample  autocorrelation  function  in  time  series  analysis.  This 
can  also  be  seen  later  in  section  3.4.5  when  the  expected  error  in  predic- 
tions is  discussed;  Ames  and  Egan's  expected  error  depends  solely  on  the 
autocorrelation  they  have  calculated.  They  also  take  this  position  when, 
after  using  a negative  exponential  function  to  approximate  the  autocorrela- 
tion, they  then  compare  it  to  their  estimate  of  p concluding  that  the 
approximation  is  only  a few  percent  larger  than  the  observed  values  of  the 
autocorrelation.  They  simply  compare  two  different  estimates  of  the  true 
state  of  nature. 

In  the  appendix  of  their  paper,  Ames  and  Egan,  [1],  give  the  following 
equation  to  compute  the  autocorrelation: 


o(t,T)  = 


lH  Cx.(t)  - x(t)]  [x.(t-T)  - x(t-T)] 

i*l  1 1 


Cx1(t)  - x(t)]2 


ZN  [x-(t-T)  - x(t-T)]2 
i=l  1 


(3.4.7) 


Of  course,  they  do  not  obtain  p(t,T)  the  true  state  of  nature,  but  just  an 
estimate  of  it,  say  r(t,T).  The  estimate  (3.4.7)  is  a function  of  the  time 
t and  the  prediction  lead  time  T. 

The  above  estimate  of  p(t,T)  is  essentially  equivalent  to  the  estimate 
given  by  the  following  expression: 


rxx(lt>  * 


t 

t«l 


N-k 


(xt  " *1  } (Vk  ‘ *2} 


C £N'k  (xt  - ZN_k  (x  . - x2)] 

t»l  Z 1 t=l  z K ‘ 


1/2 


(3.4.8) 


where  x-|  and  x„  are  the  means  of  the  first  and  the  last  (N-K)  observations, 
respectively.  Equation  (3.4.8)  expresses  the  lead  time  in  terms  of  lag  k 
where  r (k)  is  a function  only  of  the  lag  k.  The  estimates  (3.4.8)  and 
(3.4.7)  are  not  recommended  [3]  in  estimating  the  autocorrelation  of  the 
grounds  that,  although  it  gives  a reasonable  estimate  of  p when  considered 
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in  isolation  from  other  values  of  the  autocorrelation  function,  it  does  not 
give  a satisfactory  estimate  when  a set  of  estimates  is  required.  The  main 
disadvantage  of  (3.4.8),  and  hence  (3.4.7),  is  that  two  means  are  used  for 
the  mean  correction  and  that  these  change  with  lag;  in  addition,  the  normal 
izing  factor  changes  with  lag  k.  The  net  result  of  these  modifications  is 
that  these  estimates  are  not  positive  definite  which  violates  the  property 
of  the  autocorrelation  function. 

We  suggest  the  following  estimate,  [3],  [6  ],  which  gives  the  most 
satisfactory  estimate  of  the  sample  autocorrelation  function: 


xx 


cxx(k) 

(k,s^ 


(3.4.9) 


where 

cxx(k)  = jjf  ^~k  Ut  - x)  (xt+k  - x),  k * 0,  1,  ....  N-l  . 

(3.4.10) 

(cxx(k)  is  the  autocovariance  function).  Equations  (3.4.9)  and  (3.4.10)  are 
functions  only  of  the  lag  k and  are  independent  of  time. 


3.4.4  On  the  Prediction  Process 

Recall  that  Ames  and  Egan  use  the  following  equation  to  forecast  future 
value  of  the  ionospheric  data: 

x(t+T)  = x(t+T)  +>  [x(t)  - x ( t ) ] p(t,T)  . (3.4.11) 

Further,  recall  how  the  ionospheric  data  is  collected.  This  is  essentially 
what  occurs:  first  a day  is  divided  into  144-  10-minute  increments;  then 
ionospheric  soundings  are  transmitted,  then  received  and  recorded  for  each 
10-minute  increment  throughout  the  day;  this  is  repeated  for  several  days. 

The  ionospheric  series,  unlike  the  series  usually  encountered  in  practice, 
have  several  realizations  for  each  time  slot.  Initially,  to  begin  the  predic- 
tion process  using  (3.4.11),  it  is  necessary  to  speculate  or  reckon  the 
predictions  for  the  first  four  or  five  days  until  sufficient  information 
becomes  available.  Once  four  or  five  days  of  data  have  been  collected. 
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"fairly"  smooth  estimates  of  x(t+T),  x(t),  and  p(t,T)  can  be  obtained  and  the 
prediction  process  (3.4.11)  really  begins.  We  question  how  good  the  esti- 

A 

mates  of  x(t+T),  x(t),  and  p(t,T)  are  from  samples  of  only  four  or  five 
points.  In  view  of  this,  we  are  lead  to  the  statement  by  Yaglom,  [4  ],  "If 
the  number  of  observed  values  of  a series  is  small  (<  10  say),  then  the 
entire  formulation  of  the  problem  is  clearly  quite  unrealistic,  since  we 
would  not  be  able  to  make  a sufficiently  reliable  determination  of  the  auto- 
correl ati on  functi on . " 

Also,  we  feel  that  the  prediction  model  (3.4.11)  is  quite  unrealistic 
with  respect  to  the  expression  x(t+T).  It  would  appear  that  Ames  and  Egan 
use  the  value  they  are  trying  to  predict  in  calculating  the  mean  of  the  t+T 
time  slot.  The  paper  was  extremely  vague  in  how  this  mean  was  calculated. 

We  question  the  use  of  Fourier  coefficients  in  converting  and,  hence, 
in  reducing  the  data  for  computer  memory  capacity  for  two  reasons.  First, 
we  fail  to  see  how  the  required  means,  autocorrelation,  etc.,  can  be  removed 
so  easily  from  the  Fourier  coefficients  as  needed.  Secondly,  the  ionospheric 
data  appears  to  exhibit  non-stationary  properties  (this  has  been  verified  by 
the  examination  of  other  ionospheric  data),  and  Fourier  analysis  breaks  down 
when  applied  to  data  which  exhibit  random  changes  of  frequencies,  amplitudes, 
and  phases,  since  Fourier  analysis  is  based  on  the  assumption  of  fixed 
frequencies,  amplitudes,  and  phases. 

We  feel  that  the  prediction  process  of  A.,ies  and  Egan  is  quite  limited 
in  several  respects  due  in  part  to  the  criticisms  presented  above.  First, 
very  few  practical  problems  (in  time  series)  arise  such  that  we  would  have 
more  than  one  realization  for  a specific  time  slot.  Ames  and  Egan's  predic- 
tion process  (3.4.11)  is  developed  solely  from  the  viewpoint  of  having 
several  realizations  for  each  time  slot.  Secondly,  as  we  mentioned  in  the 
previous  paragraphs,  it  is  quite  unrealistic  to  estimate  the  parameters  in 
the  prediction  model  with  so  small  a sample.  (Ames  and  Egan  make  no  refer- 
ence as  to  how  much  previous  data  is  necessary  when  utilizing  their  approach). 
Finally,  it  seems  very  expensive  to  require  so  many  previous  days  of  data 
(at  least  10)  in  order  to  properly  employ  the  model. 
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3.4.5  Expected  Error  in  Predictions 

Ames  and  Egan  state  that  the  expected  error  depends  on  the  autocorrela- 
tion and  is  given  by: 

o (t+T)  « o(t+T)  / 1 - p2 (t,T)  . (3.4.12) 

/\ 

a(t+T)  is  the  sample  standard  deviation  computed  for  the  t+T  time  slot. 

Under  the  premise  that  the  expected  error  depends  on  the  autocorrelation, 
they  are  able  to  conclude  that  the  expected  reduction  of  variance  is  achieved 

A 

when  p is  greater  than  0.86. 

Then,  to  measure  the  performance  of  the  prediction  model  (3. 4. Tl), they 
used  the  following  expression: 

- HGF(t E /Wt)  (3.4., 3) 

and  compared  the  results  to  (3.4.12).  From  this  comparison,  Ames  and  Egan 
were  able  to  conclude  that  the  prediction  process  performed  nearly  as  well 
as  was  theoretically  expected.  To  have  the  theorectically  expected  error, 
it  would  have  been  necessary  to  have  the  true  states  of  nature  in  equation 
(3.4.12).  It  is  apparent  that  Ames  and  Egan  have  again  assumed  the  errone- 
ous position  that  they  have  the  true  states  of  nature,  whereas,  they  have 
only  estimates. 

It  appears  that  the  expected  error  would  depend  to  a certain  extent  on 
the  autocorrelation  function  due  to  its  presence  in  the  prediction  model 
(3.4.11).  However,  Ames  and  Egan  are  quite  vague  on  this  point  and  give  no 
derivation  of  the  dependence  of  the  expected  error  on  the  autocorrelation, 
equation  (3.4.12). 

It  seems  that  a better  criterion  to  determine  the  gordness-of-fit  of 
a fitted  model  or  a prediction  process  would  be  given  by  the  squared  error 
loss;  that  is, 

M.S.E.  » 2N  [x(t)  - x(t)]2  , (3.4.14) 

t=l 

where  x(t)  is  the  observed  series  and  x(t)  is  the  predicted  (modeled)  series. 
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Regardless  of  the  criticisms  that  we  have  made  in  this  section,  some 
interesting  ionospheric  predictions  were  obtained  using  Ames  and  Egan's 
prediction  process.  In  table  3.13,  a comparison  is  made  between  the  best 
predictions  obtained  from  the  use  of  Ames  and  Egan's  model,  made  for  the  13th 
day  of  the  experiment,  for  lead  times  of  1 , 2 and  3 steps  ahead  against  the 
time-series  approach  without  updating. 


Table  3.13  A Comparison  of  Forecasts  for  the  13th  Day  01 
Observations  at  Origin  t = 72  Between  the 
Ames-Egan  and  Time  Series  Approaches 


TIME 

OBSERVED 

VALUE 

LEAD 

TIME 

AMES-EGAN 

FORECAST 

DIFFERENCE 

TIME  SERIES 
FORECAST 

DIFFER- 

ENCE 

1150 

9.30 

1 

9.12 

-.18 

9.17 

-.13 

9.60 

2 

9.14 

-.46 

9.13 

-.47 

1210 

9.60 

3 

10.09 

+ .49 

9.17 

-.43 

The  t = 72  origin  occurs  at  the  most  stable  time  of  day  for  the  ionosphere. 
Therefore,  from  the  Ames-Egan  point  of  view,  these  forecasts  are  among  the 
best  possible  over  the  24-hour  period.  Other  times  of  day  yield  much  poorer 
forecasts  with  their  method.  Clearly,  the  time- series  approach  is  better  not 
only  from  the  theoretical  point  of  view,  as  we  have  shown,  but  also  from  the 
actual  excercising  of  the  models. 

3.5  SUMMARY  AND  CONCLUSIONS 

In  this  section,  the  procedural  approach  developed  in  Section  2.4  was 
followed  precisely  in  characterizing  actual  data.  Namely,  ionospheric  data 
obtained  by  sounding  the  ionosphere  in  the  HF  range  was  modeled  and  analyzed. 
Specifically,  the  experimental  design  and  acquisition  of  data  over  the  500  Km 
path  between  Ft.  Monmouth,  N.  J.,  and  Ft.  Drum,  N.  Y.,  was  discussed.  The 
resulting  information,  which  was  time  dependent,  consisted  of  collecting  VI 
and  01  soundings  every  ten  minutes  throughout  a 24-hour  period.  We  justified 
that  the  data  were  indeed  non-stationary  stochastic  realizations,  and  then 
proceeded  to  perform  a time  series  analysis.  The  thrust  of  this  section  was 
towards  analyzing  four  stochastic  realizations,  two  of  which  were  randomly 
selected  VI  and  01  diurnal  series,  and  two  of  which  were  the  18-day  averaged 
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VI  and  01  series.  The  random  selection  process  yielded  the  13th  day 
ionospheric  observations  for  analysis  along  with  the  averaged  information. 

Following  the  procedural  approach  recommended  in  Section  2.4,  the 
following  stochastic  processes  were  formulated  as  the  most  appropriate 
characterizations  of  the  given  information: 

i.  xt  = 1.289xt_1+  0.106xt_2  - 0.395xt_3  + .007  + Zt  for  the  13th 

day  observed  VI, 

ii.  xt  ■ 1.324xt_T  - 0.207xt_2  + 0.184xt_3  - 0.043xt_4  - 0.258xt_5 

+ .001  + for  the  18-day  averaged  VI. 

iii.  xt  = 1.484xt_1  - 0.586xt_2  + 0.368xt_3  - 0.266xt_4  + Z t for  the 

13th  day  oo served  01, 

iv.  x^  = 1.426xt_-j  - 0.204xt_2  - 0.060xt_3  - 0.212xt_4  - 0.070xt_3 

+ .0016  + Z^  f8r  the  18-day  averaged  01. 

In  selecting  these  models,  we  utilized  the  criterion  of  minimum  residual 
variance  because,  as  indicated  in  Section  2,  we  believe  this  to  be  the  most 
appropriate  criterion  for  decision  with  respect  to  identifying  the  actual 
difference  equations  which  characterize  the  ionospheric  information. 
Furthermore,  we  have  structured  tables  that  show  the  short  and  long-term 
forecasts  of  VI  and  01  soundings  along  with  their  confidence  limits.  These 
models,  in  addition  to  being  useful  for  prediction  purposes,  can  be  utilized 
in  formulating  the  theoretical  spectrum.  Such  a spectrum  would  be  extremely 
useful  in  comparing  the  smoothed  spectral  density  of  the  raw  information  with 
respect  to  identifying  the  most  useful  spectral  density  estimate  which  will 
convey  information  concerning  the  distribution  of  variance  as  a function  of 
time.  Such  information  will  be  useful  in  designing  more  efficient  HF  commun- 
ications systems. 

In  addition,  we  have  discussed  the  Ames-Egan  model  for  predicting 
ionospheric  conditions.  The  shortcomings  concerning  the  relevance  of  this 
model  were  discussed  in  some  detail  in  section  3.4. 
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4.  SPECTRAL  ANALYSIS  OF  VERTICAL  INCIDENCE 

and 

SHORT-PATH  OBLIQUE  INCIDENCE  IONOSPHERIC  SOUNDINGS 
4.1  INTRODUCTION 

In  Section  3.3,  a detailed  modeling  procedure  was  illustrated  that  yielded 
appropriate  time  series  models  for  the  ionospheric  data  described  in  section 
3.2.  It  was  concluded  that  the  models  obtained  did  characterize  the  true 
underlying  stochastic  processes  to  a high  degree.  However,  additional  anal- 
ysis of  this  information  is  necessary  to  utilize  the  ionospheric  media  more 
efficiently.  Additional  information,  therefore,  will  be  sought  with  regard 
to  the  distribution  of  the  variance  of  the  filtered  data  with  respect  to 
frequency.  Thus,  we  will  utilize  the  power  spectra  to  describe  in  detail  how 
the  variance  of  the  non-stationary  realizations  are  distributed  with 
frequency  of  occurrence  (not  the  observed  critical  frequencies). 

Ionospheric  information  is  usually  collected  at  individual  stations  as 
VI  data,  and,  for  the  benefit  of  communicators  operating  over  specified  paths, 
is  usually  translated  into  equivalent  01  information  through  the  classical 
Secant  <p  Law,  [2].  Since  there  are  an  infinite  number  of  oblique  paths  that 
can  be  utilized  by  communicators,  one  can  see  the  importance  of  converting 
the  VI  information.  Also,  consider  that  the  number  of  VI  sounder  stations 
is  limited  throughout  the  world.  This  means  that  any  such  data  acquired  at 
one  station  may  also  be  translated  into  VI  information  suitable  for  interpre- 
tation at  other  geographical  locations  within  reason.  Relationships  for 
various  translations  that  have  been  developed  in  the  United  States  by  the 
National  Bureau  of  Standards,  [2],  and  used  by  communicators  throughout  the 
world,  can  be  traced  back  to  at  least  1941. 

The  classical  Secant  <t>  Law  is  a widely  used  linear  relationship  between 
01  and  VI  data.  With  Secant  $ used  as  the  obliquity  factor,  the  law  simply 
stated  Is: 

Xqj  a Xyj  ' ^ * 

where  Xyj  is  the  observed  vertical  Incidence  information,  <p  is  the  angle  of 
incidence  of  the  radio  wave  path  at  its  entrance  into  the  ionosphere,  and 
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X0I  is  the  equivalent  information  over  an  oblique  path.  The  linear 
relationship  certainly  cannot  be  expected  to  yield  suitable  equivalent  01 
information  in  view  of  the  highly  stochastic  nature  of  the  ionosphere.  The 
Secant  <p  Law  relies  on  the  assumptions  of  a spectral  reflection  of  energy 
from  the  ionosphere  (actually,  reflections  are  dispersive  in  nature)  and  on 
the  concept  of  stratified  ionospheric  layers.  This  implies  homogeneous  iono- 
spheric layers,  where  boundaries  between  layers  of  electron  density  are  always 
definable.  This,  of  course,  is  not  the  case  since  the  ionosphere  is  random 
and  inhomogeneous,  and  it  is  affected  by  a variety  of  anomalous  activity, 
i.e. , sunspots,  magnetic  storms,  diurnal  and  seasonal  changes  in  structure. 
The  Secant  0 Law,  [2],  also  implies  that  as  path  distance  increases,  the 
oblique  information  becomes  more  uncorrelated  with  the  vertical  incidence 
information.  At  the  500  Km  path  distance,  over  the  specified  Fort  Monmouth 
- Fort  Drum  path,  the  difference  between  01  and  VI  became  significant  with 
respect  to  forecasting  ionospheric  conditions  over  the  path  using  VI  data 
alone. 

Thus,  one  can  see  that  additional  information  as  to  the  distribution  of 
the  variance  of  the  filtered  01  and  VI  information  is  extremely  important, 
and  that  information  on  the  bivariate  behavior  of  the  t o is  essential  in 
order  to  gain  a more  realistic  view  of  the  relationship  between  VI  and  01. 

In  the  succeeding  sections,  a detailed  spectral  analysis  of  the  13th  day 
observed  01  and  VI  information  will  be  performed.  In  section  4.2,  the  basic 
concept  of  "aliasing"  will  be  presented.  The  univariate  spectral  analysis 
will  be  addressed  In  section  4.3.  The  bivariate  analysis,  which  consists  of 
co-spectral,  quadrature  spectral,  cross-amplitude  spectral,  phase,  and 
coherency  spectral  estimates,  will  be  presented  in  section  4.4.  A summary 
and  conclusions  are  given  in  section  4.5. 

4.2  BASIC  CONCEPTS  OF  "ALIASING" 

With  regard  to  the  three  windows  described  in  detail  in  section  2.9,  the 
concept  of  aliasing  is  predicated  on  the  suppositions  that: 

i)  the  window  should  not  be  too  wide,  exposing  a significant  amount  of 
disturbances  (peaks  and  valleys  of  the  power  spectra),  and, 

ii)  at  the  same  time,  the  window  should  not  be  closed  too  far  so  as  to 
avoid  seeing  the  disturbances. 
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A significant  point  is  that  one  should  utilize  ingenuity  to  detect  the 
appropriate  window  length,  L,  so  that  irrelevant  information  is  not  exposed 
and  relevant  information  is  not  eliminated.  One  can,  by  visually  inspecting 
the  power  spectra  as  L is  varied,  determine  the  optimum  L that  will  play  the 
critical  role  in  interpretation  of  the  final,  smoothed  spectral  estimates. 

As  L increases,  the  associated  confidence  intervals  decrease.  Since  we 
wish  to  minimize  the  confidence  interval,  so  as  to  more  closely  approach  the 
true  state  of  nature,  a precise  analysis  is  done  for  each  window  with  respect 
to  L.  This  is  done  by  using  the  theoretical  spectrum  as  a guide  to  the 
relevant  peaks  and  valleys.  Having  analyzed  each  of  the  lag  windows,  i.e., 
those  of  Bartlett,  Tukey,  and  Parzen  (refer  to  section  2.9),  the  decision  is 
then  made  as  to  the  most  appropriate  L to  be  used  for  the  filtered  data.  As 
a result  of  our  analysis,  the  most  appropriate  L (each  window  may  have  a 
different  optimal  L)  is  the  one  that  identifies  itself  with  the  optimal 
window. 

In  the  following  sections,  we  will  utilize  this  philosophy  and  the 
detailed  procedure  of  sections  2.5  through  2.9  to  obtain  the  univariate  and 
bivariate  spectral  estimates. 

4.3  UNIVARIATE  SPECTRAL  ANALYSIS  OF  IONOSPHERIC  INFORMATION 

In  this  section,  an  analysis  of  two  univariate  time  series  will  be 
presented.  Specifically,  the  series  corresponding  to  the  13th  day  VI  and  01 
information  (refer  to  section  3.2)  will  be  analyzed  by  the  method  of  power 
spectra,  described  in  sections  2.5  through  2.9.  The  ionospheric  data  was 
obtained  every  10  minutes  throughout  the  day  (24  hours)  for  a total  of  144 
recordings.  Estimates  of  the  spectral  density  function  are  obtained  from  the 
filtered  data  using  the  Bartlett,  Tukey,  and  Parzen  lag  windows  described  in 
section  2.9. 

4.3.1  Estimate  of  the  Spectral  Density  Function  Using  Bartlett's  Lag  Window 

The  values  of  the  estimate  of  the  spectral  density  function  using 
Bartlett's  lag  window,  equation  (2.9.4),  were  calculated  and  plotted  versus 
frequency  for  L * 8,  12,  16,  24,  and  32  units.  As  a basis  of  comparison,  the 
estimates  were  plotted  on  the  same  set  of  axes.  For  these  values  of  L,  and 
for  L * 20,  we  calculated  the  bandwidth,  the  confidence  intervals,  and  the 
degrees  of  freedom  which  are  shown  in  table  4.1. 
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Table  4.1  Truncation  Point,  Bandwidth,  Degrees  of  Freedom,  and 
Confidence  Intervals  for  Bartlett’s  Lag  Window 


L 

Bandwidth 

d.f . 

95%  C.I. 

for  Log  Txx(f) 

8 

.188 

54 

-.154 

.176 

12 

.125 

36 

-.179 

.226 

16 

.094 

27 

-.204 

.267 

20 

.075 

22 

-.223 

.301 

24 

.063 

18 

-.243 

.340 

32 

.047 

13 

-.279 

.414 

The  formula  used  for  the  bandwidth  of  the  estimate  of  the  spectral  density 
function  is  given  by: 


1 _ 1.5 
'LA  ~ T" 


and  the  equation  for  the  degrees  of  freedom  is  given  by: 

v = 2 £ b-,=  2Tb  =»  2 ( 1 44 ) b = 288b  . 

Note  that  since  we  have  chosen  A = 1,  we  have  L * M. 

Figures  4.1  and  4.2  give  a comparison  of  the  theoretical  spectral 
density  functions  of  the  AR  processes  which  characterize  the  VI  and  01  series, 
respectively  (equations  3.3.1  and  3.3.3),  and  the  smoothed  estimates  for  the 
various  truncation  points,  along  with  the  95%  confidence  intervals. 

It  is  a known  fact  that  increasing  the  bandwidth  of  the  estimate  of  the 
spectral  density  involves  increasing  the  amount  of  bias  and  decreasing  the 
variance;  thus,  a compromise  has  to  be  reached  as  to  which  is  the  best  value 
of  L.  In  making  such  a decision,  one  should  take  into  consideration  the 
confidence  interval,  the  degrees  of  freedom,  and  the  visual  appearance  of  the 
plot  of  the  estimates.  For  L = 8,  the  plots  are  very  smooth  and  have  a shape 
which  follows  the  trend  of  the  theoretical  spectrum  with  the  bandwidth  being 
wide  enough  to  conceal  any  peaks  that  may  be  present.  In  both  cases,  by  increas- 
ing L to  12,  we  obtain  an  indication  of  other  peaks  appearing  at  f * .20  and 
.44  cycles  per  second  for  the  01  spectrum  and  at  f * .20,  .30,  and  .39  cycles 
per  second  for  the  VI  spectrum.  These  are  in  addition  to  the  major  peaks  in  the 
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FIGURE  4.1  ESTIMATE  OF  THE  SPECTRAL  DENSITY  OF  THE  FILTERED 
VI  DATA  USING  THE  BARTLETT  LAG  WINDOW 
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ESTIMATE  OF  THE  SPECTRAL  DENSITY  OF  THE  FILTERED 
01  DATA  USING  THE  BARTLETT  LAG  WINDOW 
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theoretical  spectral  densities.  The  plots  are  still  quite  smooth  and  the  j 

bandwidth  is  wide  enough  to  give  a great  deal  of  faith  to  the  estimates.  By 
increasing  L to  16,  the  bandwidth  seems  to  be  in  a very  shakey  range.  In  s 

this  case,  both  spectra  display  the  peaks  of  the  theoretical  spectral  densi- 
ties and  also  those  of  the  L 3 12  spectra.  However,  the  curves  have  changed 
very  little  from  those  for  which  L 3 12.  Since  larger  values  of  L produce 
many  small  erratic  peaks,  we  chose  L 3 16  to  estimate  the  spectral  densities 
of  the  VI  and  01  filtered  data  using  Bartlett's  lag  window. 

4.3.2  Estimate  of  the  Spectral  Density  Function  Using  Tukey's  Lag  Window 

" " " 

Using  Tukey's  lag  window  given  by  equation  (2.9.6)  the  smooth  spectral 
density  estimates  IT  (f)  for  the  VI  and  01  filtered  data  were  calculated  for 
L 3 8,  12,  14,  16,  and  24  units.  Figures  4.3  and  4.4  show  the  spectral 
density  estimates  and  the  theoretical  spectra  for  the  filtered  data  using  the 
Tukey  lag  window  for  the  various  truncation  lengths  along  with  the  95%  confi- 
dence intervals  and  the  various  bandwidths  associated  with  these  truncation 
points.  It  is  clear  that  for  L 3 8,  the  sample  spectra  have  the  same  general 
shape  as  the  theoretical  spectra  and  the  curves  are  very  smooth.  By 
increasing  the  truncation  value  to  12,  the  plots  are  still  fai/ly  smooth,  but 
peaks  appear  in  both  estimates  at  about  f 3 .19,  .31,  and  .39  cycles  per 
second,  and  at  f = .21  and  .45  cycles  per  second,  respectively.  At  L = 16, 
the  peaks  are  slightly  more  pronounced,  and  as  L is  increased  above  16,  more 
peaks  appear  at  higher  frequencies.  This  indicates  that  the  variances  are 
increasing  and,  thus,  the  sample  spectra  are  becoming  more  erratic  for 
L > 16.  On  this  basis,  it  was  decided  to  try  to  obtain  better  estimates  than 
those  calculated  by  computing  additional  spectral  density  estimates  for 
L 3 14  units. 

Table  4.2  shows,  for  the  various  truncation  points,  the  bandwidth, 
degrees  of  freedom,  and  confidence  intervals  using  Tukey's  lag  window.  The 
bandwidth  b 3 1.33/L. 
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FIGURE  4.3  ESTIMATE  OF  THE  SPECTRAL  DENSITY  OF  THE  FILTERED 

VI  DATA  USING  THE  TUKEY  LAG  WINDOW 
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RE  4.4  ESTIMATE  OF  THE  SPECTRAL  DENSITY  OF  THE  FILTERED 
01  DATA  USING  THE  TUKEY  LAG  WINDOW 


Table  4.2  Truncation  Point,  Bandwidth,  Degrees  of  Freedom,  and 
Confidence  Intervals  for  Tukey's  Lag  Window 


L 

Bandwidth 

d.f. 

95%  C.I. 

for  Loq  Txx(f) 

8 

.166 

47 

-.159 

.195 

12 

.111 

31 

-.192 

.246 

14- 

.095 

27 

-.204 

.267 

16 

.083 

23 

-.219 

.294 

20 

.067 

19 

-.238 

.329 

24 

.055 

15 

-.263 

.380 

32 

.042 

12 

CO 

CO 

C\l 

r 

.436 

Table  4.2  is  quite  helpful  in  deciding  that  for  L s 14  units,  we  will  have 
the  best  estimate  of  the  spectrum  using  Tukey's  lag  window.  The  degrees  of 
freedom,  v = 27,  are  sufficient  for  fairly  small  95%  confidence  intervals, 
and  this  gave  a bandwidth  of  .095  so  that  peaks  in  the  time  spectrum  of  band- 
widths  larger  than  .095  will  be  detected.  Decreasing  the  bandwidth  to  .083, 
that  is,  L = 16,  causes  a loss  of  four  degrees  of  freeom  and  a slight 
increase  in  the  confidence  interval  width.  For  L * 12,  the  bandwidth  is 
considerably  larger  (.111),  and  there  is  not  much  change  in  the  confidence 
interval  even  though  there  are  31  degrees  of  freedom.  Therefore,  for  a 
truncation  length  of  14  units,  we  obtain  the  best  estimate  for  both  spectra 
using  Tukey's  lag  window. 

As  we  mentioned  previously,  the  plot  of  the  estimate  of  the  spectral 
density  is  given  ir.  the  logarithmic  scale  to  show  more  detail  in  the  spectrum 
over  a wider  amplitude  range. 

4.3.3  Estimates  of  the  Spectral  Density  Function  Using  Parzen's  Lag  Window 

Using  Parzen's  lag  window,  given  by  equation  (2.9.8),  we  obtained  esti- 
mates of  the  spectral  density  functions  for  the  two  first  order  filtered  VI 
and  01  data,  for  various  truncation  points.  As  before,  we  shall  let  A * 1, 
so  that  L = M,  the  truncation  points  of  the  smoothed  spectral  estimator.  L 
was  varied  from  8 to  32  in  intervals  of  four  and  eight  units. 

Figures  4.5  and  4.6  show  the  spectral  density  estimates  of  both  the  VI 
and  01  filtered  series  for  the  various  truncation  points  along  with  the 
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FIGURE  4.5  ESTIMATE  OF  THE  SPECTRAL  DENSITY  OF  THE  FILTERED 

VI  DATA  USING  THE  PARZEN  LAG  WINDOW 
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FIGURE  4.6  ESTIMATE  OF  THE  SPECTRAL  DENSITY  OF  THE  FILTERED 

01  DATA  USING  THE  PARZEN  LAG  WINDOW 


theoretical  spectral  density  of  the  respective  autoregressive  processes. 

In  addition,  95%  confidence  intervals  and  the  corresponding  bandwidths  are 
displayed.  The  bandwidth  using  Parzen's  lag  window  is  given  by: 

b = 1^§6  . K86 

The  degrees  of  freedom  for  the  confidence  intervals  were  found  using  the 
following  relationship: 


V = 2 J b1  * 288b  , 

where  b-j  = 1.86  for  the  Parzen  window  and  T = total  of  observations,  which 
in  these  cases  is  144  soundings.  Table  4.3  gives,  for  the  various  truncation 
points,  the  corresponding  bandwidths,  degrees  of  freedom,  and  a 95%  confi- 
dence interval  for  the  theoretical  spectrum,  Txx(f),  for  the  Parzen  lag 
window. 


Table  4.3 

8andwidth, 
for  Selected 

Degrees  of  Freedom,  and 
Values  of  L for  Parzen' 

95%  Confidence  Intervals 
s Window 

L 

Bandwidth 

d.f . 

95%  C.I. 

for  Loq  Txx(f) 

8 

.233 

67 

-.135 

.160 

12 

.155 

44 

-.164 

.203 

16 

.116 

33 

-.186 

.235 

20 

.093 

27 

-.204 

.267 

24 

.078 

22 

-.223 

.301 

32 

.058 

16 

-.255 

.365 

In  selecting  a proper  value  for  L for  the  spectral  densities,  one  should 
be  able  to  detect  peaks  in  the  spectra,  have  reasonable  confidence  intervals, 
and  a bandwidth  that  provides  a reasonable  bias.  For  an  l value  of  8 units, 
the  spectral  densities  were  much  too  smooth,  and  we  were  unable  to  detect 
peaks  less  than  .233  wide.  Increasing  the  L values  from  16  to  20  units  gave 
a fairly  reasonable  display  of  both  spectral  densities.  At  L * 20,  two  major 
peaks  occur  that  are  quite  similar  to  those  of  the  respective  theoretical 
densities.  For  a truncation  point  of  24  units,  very  small  peaks  began  to 
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appear  which  indicate  that  the  variances  may  be  influencing  each  of  the 
densities.  This  was  also  evident  at  L = 32,  where  the  peaks  became  very 
erratic  and  very  noticeable.  Thus,  the  choice  was  narrowed  very  quickly  to 
choosing  either  L = 16  or  20  units.  The  confidence  intervals  for  L s 16  and 
20  are  .421  and  .471,  respectively.  The  bandwidth  for  L * 20,  however,  is 
smaller  by  about  20%  from  that  of  L 5 16.  Therefore,  the  spectral  densities 
corresponding  to  L = 20  units  were  selected  as  the  most  reasonable.  The 
spectral  density  estimates  clearly  show,  for  both  the  VI  and  01  filtered 
series,  that  most  of  the  power  is  concentrated  at  low  frequencies.  For 
example,  a major  peak  for  the  01  spectrum  is  located  at  f * .15  cycles  per 
second  with  smaller  peaks  located  at  f 5 .30  and  .41  cycles  per  second.  The 
bandwidth  for  L = 20  units  is  .093,  which  means  that  we  can  detect  peaks  with 
a width  of  this  value  or  greater.  The  above  remarks  are  graphically  verified 
in  figures  4.5  and  4.6  where  the  theoretical  spectral  density  for  the  respec- 
tive autoregressive  models  are  compared  with  the  spectral  estimate  for  L = 3, 
12,  16,  20,  and  24  units.  In  addition,  the  95%  confidence  intervals  and  the 
corresponding  bandwidths  for  the  truncation  points  are  given. 

4,4  BIVARIATE  SPECTRAL  ANALYSIS  OF  THE  IONOSPHERIC  INFORMATION 

In  this  section,  we  shall  be  concerned  with  analyzing  the  bivariate 
behavior  of  the  13th  day  observed  vertical  and  oblique  incidence  ionospheric 
soundings  for  the  500  Km  experiment.  More  specifically,  estimates  of  the 
smoothed  quadrature,  phase,  and  cross-amplitude  spectra  will  be  obtained 
using  the  three  lag  windows  discussed  in  section  2.9.  In  addition,  estimates 
of  the  coherency  spectrum  will  be  obtained. 

Having  calculated  and  plotted  the  cross-amplitude  spectrum,  one  can 
detect  whether  or  not  frequency  components  in  the  vertical  incidence  sound- 
ings are  associated  with  large  or  small  amplitudes  at  the  same  frequency  in 
the  01  series.  The  estimate  of  the  phase  spectrum  of  the  two  stochastic 
realizations  helps  us  in  determining  whether  or  not  frequency  components  in 
the  VI  series  are  in  phase  or  out  of  phase  (lag  or  lead)  with  components  at 
the  same  frequency  in  the  01  series.  The  cross-amplitude  spectrum,  Ayz(f), 
is  a measure  of  the  covariance  that  exists  between  the  01  and  VI  soundings 
at  frequency  f,  and  is  the  square  root  of  the  sum  of  the  squares  of  the 
co-spectral  and  coquadrature  spectral  estimates.  An  estimate  of  the 
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cross-amplitude  spectrum  and  the  'phase  spectrum  would  suffice  to  provide  a 
complete  description  of  the  behavior  of  the  two  series. 

In  general,  the  squared  coherency  spectrum  plays  the  role  of  a correla- 
tion coefficient  with  respect  to  frequency.  Its  usefulness  lies  in  the  fact 
that  dimensions  do  not  enter  the  picture  when  the  correlation  is  measured 
with  respect  to  frequency.  Unlike  the  squared  coherency  spectrum,  the  cross- 
amplitude spectrum  depends  upon  the  dimensions  of  the  01  and  VI  soundings. 

This  is  th«»  reason  the  squared  coherency  spectrum  is  sometimes  preferred  over 
the  cross-amplitude  spectrum  and,  together  with  the  phase  spectrum,  will  give 
a complete  picture  of  the  cross  correlation  behavior  of  the  01  and  VI 
soundings. 

With  respect  to  the  aims  of  the  present  study,  we  will  only  give  the 
equations  (estimates)  that  characterize  the  above  concepts,  and  we  will  not 
discuss  the  theoretical  implications.  For  complete  details  of  these  concepts, 
refer  to  sections  2.5  through  2.9. 

4.4.1  Co-Spectrum  Estimates  Using  the  Bartlett,  Tukey,  and  Parzen  Lag 
Windows 

We  shall,  in  what  follows,  obtain  estimates  for  the  co-spectral , quadra- 
ture, phase,  and  cross-amplitude  spectral  estimates  using  Bartlett’s  lag 
window.  These  smoothed  estimates  were  obtained  using  the  truncation  points 
L 3 M 3 3,  12,  16,  24,  and  32  units.  These  truncation  points  correspond  to 
decreasing  the  bandwidth  to  b 3 b^/L  = 1.5/L. 

Figure  4.7  shows  the  smoothed  co-spectral  estimates.  Similarly,  figure 
4.8  shows  the  various  smoothed  quadrature  spectral  estimates.  It  is  clear 
that  for  L > 24  units,  the  estimates  in  both  cases  become  very  erratic.  As 
we  mentioned  previously,  compromising  between  bias  and  variance,  it  appears 
. that  for  L = 16  units,  we  have  the  best  estimate  using  Bartlett's  lag  window 
with  b 3 .094  and  v 3 27  degrees  of  freedom.  The  smoothed  cross-amplitude 
. spectral  estimate  and  the  smoothed  phase  spectral  estimate,  plotted  for 

L 3 16,  each  on  separate  sets  of  axes  to  enhance  the  details  of  the  series, 
are  shown  in  figures  4.9  and  4.10,  respectively. 

The  smoothed  co-spectral,  quadrature,  phase,  and  cross-amplitude  spec- 
tral estimates  were  similarly  obtained  using  Tukey' s lag  window.  The 


truncation  points  used  for  the  co-spectral  and  quadrature  spectral  estimates 
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FIGURE  4.7  SMOOTHED  CO-SPECTRAL  ESTIMATES 
USING  THE  BARTLETT  LAG  WINDOW 


FIGURE  4.8  SMOOTHED  QUADRATURE  SPECTRAL  ESTIMATES 
USING  THE  BARTLETT  LAG  WINDOW 
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FIGURE  4.9  SMOOTHED  CROSS- AMPLITUDE  SPECTRAL  ESTIMATE  USING  THE 

BARTLETT  LAG  WINDOW  FOR  L = 16 
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FIGURE  4.10  SMOOTHED  PHASE  SPECTRAL  ESTIMATE  USING 
THE  BARTLETT  LAG  WINDOW  FOR  L = 16 


were  L 3 8,  12,  14,  16,  20  and  24  units.  Figure  4.11  displays  the  smoothed 
co-spectral  estimates.  The  smoothed  quadrature  spectral  estimates  are 
plotted  in  figure  4.12  for  the  same  truncation  points.  For  both  of  these 
cases,  the  estimates  become  more  erratic  as  L is  increased  beyond  20  units. 
Taking  the  bandwidth  into  consideration,  we  choose  the  estimate  for  which 
L 3 14  units  as  the  best  compromise  between  bias  and  variance.  Thus,  for 
L 3 14,  the  bandwidth  resulted  in  b = 1.33/L  3 .095  and  v = 27  degrees  of 
freedom  for  the  Tukey  lag  window.  Decreasing  b to  .083,  the  degrees  of 
freedom  are  decreased  considerably,  therefore,  L 3 14  units  will  give  the 
best  estimate  of  the  co-  and  quadrature  spectra  for  the  Tukey  lag  window. 

The  smoothed  phase  and  smoothed  cross-amplitude  spectra  were  then  plotted  for 
L 3 14  units  to  enhance  the  details.  Figures  4.13  and  4.14  display  the 
smoothed  cross-amplitude  spectral  estimate  and  the  smoothed  phase  spectral 
estimate,  respectively,  using  the  Tukey  lag  window  for  L 3 14. 

A similar  analysis  was  performed  to  obtain  smoothed  estimates  for  the 
co-  and  quadrature  spectra  using  Parzen 's  lag  window  for  L = 8,  12,  16,  20 
ana  24  units.  Figures  4.15  and  4.16  display  the  above  smoothed  estimates. 

The  bandwidths  for  the  Parzen  lag  window  are  given  by  b 3 1.86/L  and  the 
degrees  of  freedom  can  be  obtained  from  v = 288b.  For  values  of  L _>  24,  the 
estimates  become  somewhat  erratic  and  the  bandwidth  and  degrees  of  freedom 
are  decreased.  However,  the  decrease  in  bandwidth  from  .093  to  .078 
(co-spectral  estimate)  for  L 3 20  and  24  units,  respectively,  is  not  worth 
the  decrease  in  variance.  Hence,  we  choose  L 3 20  as  our  best  estimates  of 
the  co-  and  quadrature  spectra.  This  gives  a bandwidth  of  b 3 .093.  Figures 
4.17  and  4.18  show  the  smoothed  cross-amplitude  and  phase  spectral  estimates, 
respectively,  for  l = 20,  using  the  Parzen  lag  window,  along  with  the  corres- 
ponding bandwidth.  Note  that  since  F-^f)  1S  sma^  » L^(f)  ~ A-^f). 

4.4.2  Choosing  the  Best  Lag  Window  and  Truncation  Point 

To  compare  the  estimates  obtained  for  the  Bartlett,  Tukey,  and  Parzen 
lag  windows,  the  estimates  corresponding  to  the  best  value  of  L (chosen  for 
each  window)  were  plotted  on  the  same  axes  (see  figures  4.19  and  4.20).  The 
estimates  for  the  co-spectra  coincided  almost  exactly.  Each  estimate  has 
27  degrees  of  freedom  for  the  autospectrum  analysis.  The  Parzen  lag  window 
has  a slightly  smaller  bandwidth  than  the  others.  It  was  difficult  to  choose 
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FIGURE  4.11  SMOOTHED  CO-SPECTRAL  ESTIMATES 
USING  THE  TUKEY  LAG  WINDOW 


FIGURE  4.12  SMOOTHED  QUADRATURE  SPECTRAL  ESTIMATES 

USING  THE  TUKEY  LAG  WINDOW 
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riGURE  4.13  SMOOTHED  CROSS- AMPL ITUDE  SPECTRAL  ESTIMATE  USING 

THE  TUKEY  LAG  WINDOW  FOR  L = 14 
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SMOOTHED  CO-SPECTRAL  ESTIMATES 
USING  THE  PARZEN  LAG  WINDOW 
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FIGURE  4.16  SMOOTHED  QUADRATURE  SPECTRAL  ESTIMATES 
USING  THE  PARZEN  LAG  WINDOW 
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FIGURE  4.17  SMOOTHED  CROSS- AMPLITUDE  SPECTRAL  ESTIMATE  USING 

THE  PARZEN  LAG  WINDOW  FOR  L = 20 
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FIGURE  4.18  SMOOTHED  PHASE  SPECTRAL  ESTIMATE  USING 
THE  PARZEN  LAG  WINDOW  FOR  L = 20 
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FIGURE  4.19  COMPARISON  OF  THE  SMOOTHED  CO-SPECTRAL  ESTIMATE  OF  THE 
BARTLETT,  TUKEY,  AND  PARZEN  LAG  WINDOWS 
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COMPARISON  OF  THE  SMOOTHED  QUADRATURE  SPECTRAL  ESTIMATE 
OF  THE  BARTLETT*  TUKEY,  AND  PARZEN  LAG  WINDOWS 


the  best  window,  but  since  Parzen's  lag  window  for  L = 20  units  gave  a band- 
width of  .093,  it  was  chosen  as  the  best  smoothed  estimate  of  the  co-  and 
quadrature  spectra.  The  smoothed  estimates  for  the  phase  and  cross-amplitude 
spectra  are  also  best  represented  by  the  lag  window  of  L = 20  units.  The 
smoothed  sample  co-spectral  estimate  estimates  the  covariance  due  to  the 
in-phase  components.  There  are  peaks  at  about  .20  and  .31  cycles  per  second 
which  correspond  to  the  peaks  in  the  autospectra  due  to  the  fact  that  the 
variance  is  a special  case  of  the  covariance  (see  figure  4.21).  At  frequen- 
cies less  than  .15  cycles  per  second,  the  covariance  between  the  VI  and  01 
realizations  is  fairly  large  and  constant  over  the  frequency  range  0 to  .15 
cycles  per  second.  The  variance  at  most  frequencies  in  the  autospectra  is 
fairly  large.  However,  the  covariance  distribution  of  the  in-phase  compo- 
nents of  the  filtered  ionospheric  series  is  small,  and  therefore,  the  series 
in-phase  components  are  not  very,  dependent.  The  larger  value  of  the  sample 
cospectrum  is  near  0 cycles  per  second  corresponding  to  variance  values  of 
autospectra  of  about  10  at  the  same  frequency  for  the  Parzen  lag  window. 
However,  L 3 20  units,  and  hence,  the  correlation  is  small  as  will  be  veri- 
fied by  the  squared  coherency  spectral  estimate. 

The  smoothed  quadrature  spectral  estimate  estimates  the  covariance  of 
the  out-of-phase  components  of  the  two  filtered  time  series.  This  also  shows 
that  there  is  small  covariance  between  the  out-of-phase  components  of  the  two 
filtered  series  and,  hence,  that  they  are  not  very  correlated.  The  largest 
value  of  the  estimate  is  .012  for  the  chosen  lag  window  (Parzen,  L = 20),  and 
the  smallest  value  is  -.011  (see  figure  4.22).  One  can  conclude,  therefore, 
that  there  is  little  covariance  exhibited  throughout  the  range  0 to  .50  cps., 
but  the  out-of-phase  components  vary  in  a sinusoidal  manner  at  all  frequen- 
cies. 

I 

. The  smoothed  phase  spectral  estimate  estimates  the  phase  angle  in 

radians  by  which  the  VI  filtered  time  series  leads  or  lags  the  filtered  01 
series  (see  figure  4.18).  At  frequencies  0 to  .05  cps.,  the  phases  are 
approximately  the  same  (phase  spectral  estimate  near  0).  At  frequencies 
between  .05  cps.  and  .20  cps.,  the  in-phase  components  of  the  two  time  series 
lag  the  out-of-phase  components  very  slightly.  From  .20  cps.  to  approxi- 
mately .27  cps.,  the  out-of-phase  components  lag  the  in-phase  components. 

From  .27  cps.  to  .37  cps.,  the  in-phase  is  lagging,  and  from  .37  cps.  to 
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SMOOTHED  QUADRATURE  SPECTRAL  ESTIMATE  USING 
THE  PARZEN  LAG  WINDOW  FOR  L = 20 


.5  cps.,  the  out-of-phase  components  lag  the  in-phase  components  of  the  two 
time  series.  Since  the  phases  alternate  leading,  there  is  no  reason  to  , 
assume  or  conclude  that  one  time  series  leads  or  lags  the  other  at  all 
frequencies. 

The  smoothed  cross-amplitude  spectral  estimate  shows  whether  or  not  the 
amplitude  of  the  components  at  a particular  frequency  in  one  time  series  is 
associated  with  a large  or  small  amplitude  of  the  same  order  at  the  same 
frequency  in  the  other  time  series.  The  spectral  density  of  the  autospectra 
shows  that  the  variance  is  about  10  in  both  filtered  series  so  that,  at 
frequencies  from  0 cps.  to  .15  cps.,  the  amplitude  of  the  components  of  one 
time  series  is  associated  with  corresponding  large  or  small  amplitudes  at  the 
same  frequency  in  the  other.  Again,  this  seems  to  indicate  that  the  covari- 
ance between  the  component  amplitudes  is  near  zero  at  other  frequencies. 

4.5  SUMMARY  AND  CONCLUSIONS 

A plot  (see  figure  4.19)  is  given  for  the  selected  best  estimates  of  the 
co-spectral  densities  for  each  of  the  three  lag  windows,  namely,  those  of 
Bartlett,  Tukey,  and  Parzen.  Although  the  truncation  is  different  for  each 
lag  window,  the  bandwidth,  degrees  of  freedom,  and  confidence  intervals  are 
almost  identical.  Thus,  it  is  quite  difficult  to  choose  which  lag  window 
gives  the  best  smoothed  estimate  of  the  spectral  density  function.  However, 
calculating  the  approximate  bias  for  each  of  the  above  lag  windows,  it  was 
found  that  the  bias  for  Parzen' s lag  window  is  somewhat  smaller  than  those 
of  the  Tukey  and  Bartlett  lag  windows.  That  is: 

Bp(f)  = W r2yy(f)  • 

Furthermore,  the  variance  ratio,  that  is,  the  proportional  reduction  in  vari- 
ance as  the  result  of  using  the  smoothed  estimator  as  compared  to  the  sample 
spectrum  estimate,  is  approximately- equal  to  .128.  On  the  basis  of  these 
two  criteria,  the  best  estimate  of  the  spectral  density  was  chosen  using 
Parzen's  lag  window.  In  addition,  the  bandwidth  of  this  lag  window  is 
slightly  smaller  than  those  of  the  Tukey  and  Bartlett  lag  windows.  There- 
fore, the  best  estimate  of  the  spectral  density  of  both  the  observed  VI  and 
01  soundings  was  obtained  using  Parzen's  lag  window  for  L * 20  units.  This 
value  of  L resulted  in  a 95*  confidence  Interval  width  of  .471  with  27  degrees 
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of  freedom,  and  a bandwidth  of  b = .093.  The  bandwidth  is  less  than  1/5  of 
the  total  frequency  range  over  which  the  spectral  density  function  is  esti- 
mated. Since  we  are  detecting  peaks  with  widths  of  .093  or  more,  the  peaks 
appearing  in  the  estimated  spectral  density  of  the  01  spectrum  at  frequencies 
f 3 .15,  .30,  and  .40  cps.,  are  valid  peaks  and  they  should  be  taken  into 
consideration  in  interpreting  the  behavior  of  the  observed  01  soundings  (see 
figure  4.6).  The  process  generating  the  soundings  exhibits  a large  variance 
around  these  three  frequencies  for  the  filtered  data.  These  frequencies 
approximately  coincide  with  the  most  critical  times  of  day  for  ionospheric 
support  of  HF  communications.  They  are:  the  pre-dawn  dip  (ionosphere 
stratifies  into  D,  E,  F-|  and  F2  layers  due  to  the  sun's  energy);  mid-day 
(where  stability  is  a function  of  various  natural  and  man-made  anomalies); 
and  twilight  (where  the  ionosphere  recombines  into  one  F layer).  Such  infor- 
mation should  be  taken  into  account  in  the  design  and  operation  of  an  HF 
communications  system.  Frequencies  above  f = .200  cps.  on  the  spectral 
estimates  gives  the  lowest  power,  that  is,  the  least  variance. 

The  Parzen  lag  window  for  L = 20  units  and  b = .093  was  used  to  obtain 
smoothed  estimates  of  the  co-  and  quadrature  spectra.  The  smoothed  estimates 
of  the  phase  and  cross-amplitude  spectra  were  also  obtained  using  the  same 
lag  window  and  L = 20  units. 

The  smoothed  sample  co-spectral  estimate  estimates  the  covariance  due 
to  the  in-phase  components.  There  is  a peak  at  about  .20  cps.  and  one  at 
.31  cps.  which  correspond  to  the  peaks  in  the  autospectra.  At  frequencies 
above  .150  cps.,  the  covariance  is  reasonably  small  and  approximately 
constant  over  the  frequency  range  of  .15  to  .50  cps.  The  variance  at  most 
frequencies  in  the  autospectra  is  fairly  large.  However,  the  covariance 
distribution  of  the  in-phase  components  of  the  two  filtered  series  is  small 
and  has,  due  to  the  ionospheric  series,  in-phase  components  that  are  not  very 
dependent.  The  larger  value  of  the  sample  spectra  is  near  0 cps.,  corres- 
ponding to  variance  values  of  the  autospectra  of  about  10,  at  the  same 
frequency,  using  the  Parzen  lag  window  for  L = 20  units.  Hence,  the  correla- 
tion between  the  01  and  VI  soundings  is  small  as  was  verified  by  the  squared 
coherency  spectral  estimate  (figure  4.23). 

The  smoothed  quadrature  spectral  estimate  estimates  the  covariance  of 
the  out-of-phase  components  of  the  filtered  VI  and  01  soundings.  It  showed 
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FIGURE  4.23  SMOOTHED  SQUARED  COHERENCY  FOR  THE 
PARZEN  LAG  WINDOW  FOR  L » 20 


that  the  covariance  between  the  out-of -phase  components  of  the  two  filtered 
series  is  small,  and  hence,  that  they  are  not  very  correlated.  The  largest 
value  is  .012  for  the  chosen  lag  window,  and  the  smallest  value  is  -.011. 

There  is  little  or  no  covariance  exhibited  throughout  the  range  from  0 to 
.50  cps.,  but  the  out-of-phase  components  vary  in  a sinusoidal  manner  at  all 
frequencies. 

The  smoothed  phase  spectral  estimate  estimates  the  phase  angle  in  radi- 
ans by  which  one  filtered  time  series  leads  or  lags  another.  At  frequencies 
0 to  .05  cps.,  the  phases  are  approximately  the  same;  that  is,  the  phase 
spectral  estimate  is  near  zero.  At  frequencies  between  .05  cps.  and  .20  cps., 
the  in-phase  components  of  the  two  time  series  lag  the  out-of-phase  compo- 
nents very  slightly.  From  .20  cps.  to  approximately  .27  cps.,  the  out-of- 
phase  components  lag  the  in-phase  components.  From  .27  cps.  to  .37  cps.,  the 
in-phase  components  are  lagging,  and  from  .37  cps.  to  .50  cps.,  the  out-of- 
phase components  lag  the  in-phase  components  of  the  two  ionospheric  time 
series.  Since  the  phase  is  alternately  leading,  there  is  no  reason  to  assume 
or  conclude  that  one  time  series  leads  or  lags  the  other  at  all  frequencies. 

The  smoothed  cross-amplitude  spectral  estimate  shows  whether  or  not  the 
amplitude  of  the  components  at  a particular  frequency  in  one  time  series  is 
associated  with  a large  or  small  amplitude  of  the  same  order  at  the  same 
frequency  in  another  time  series.  The  spectral  density  of  the  autospectra 
shows  that  the  variance  is  about  10  in  both  the  filtered  VI  and  01  soundings, 
so  that,  at  frequencies  from  0 cps.  to  .15  cps.,  the  amplitude  of  the  compo- 
nents of  one  time  series  is  associated  with  corresponding  large  or  small 
amplitudes  (at  the  same  frequency)  of  the  other.  Again,  this  indicates  that 
the  covariance  between  the  component  amplitudes  is  near  zero  at  other 
frequencies. 

In  order  to  obtain  a better  representation  of  the  important  peaks  and 
a confidence  interval,  the  squared  coherency  was  calculated  and  plotted  (see 
figure  4.23). 

A 95%  confidence  interval  was  obtained  using  the  following  expression: 
yyx(f)  - ±1.96  vOTpr 

* ±1.96  /20/ZTT.86')144  = ±.379  , 
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and  is  shown  on  the  graph  of  the  smoothed  squared  coherency  spectrum.  This 
squared  coherency  spectral  estimate  gives  the  correlation  between  the 
observed  VI  soundings  and  the  observed  01  soundings  for  the  500  Km  experiment 
At  low  frequencies,  we  have  very  little  correlation  between  the  two  filtered 
series  (maximum  of  ^ .004)  as  shown  by  the  expanded  scale.  One  can  conclude 
that  there  is  virtually  no  correlation.  The  fact  that,  at  all  frequencies, 
the  squared  coherency  approaches  zero  indicates  that  the  noise  level  is  high 
in  the  filtered  series  for  all  components  at  all  frequencies.  This  indicates 
high  variance  in  the  autospectra  for  the  corresponding  frequencies.  There- 
fore, we  conclude  that  the  01  and  VI  filtered  series  are  not  correlated 
within  the  0 to  .50  cps.  range. 
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